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Chapter 1 
Introduction 



Since the first realization of Bose-Einstein condensation (BEC) in trapped atomic vapors 
in 1995 [19], these systems have received increasing attention from both experimental and 
theoretical efforts. Most theoretical studies of these many-boson systems are based on the 
so-called mean-field methods which accurately describe much of the dilute condensate ener- 
getics. Systems are termed as dilute when the average interparticle distance is much larger 
than the range of the interaction. The main parameter characterizing the interaction in 
the dilute regime is the s-wave scattering length, a, and the diluteness condition can be 
expressed in terms of the density n as n\a\^ <C 1. 

Prom the viewpoint of quantum many-body physics, the trapped atomic vapors are 
somewhat peculiar. Well above the critical point of condensation, the gases are extremely 
dilute, and their description as non-interacting bosons is very accurate. As the condensation 
sets in, the trapped atoms are strongly compressed in real space. This makes it much 
more likely that the individual particles are within interaction range of each other, and 
interactions suddenly become very important. As a consequence the motion of the particles 
becomes correlated and both the order of these correlations, that is, the number of particles 
which are simultaneously within interaction range, and their overall infiuence will depend 
on the size of n\af. 

For the density ranges attained in BEC experiments, the diluteness condition may well 
be broken, exploiting the large variation of the scattering length in the vicinity of a Feshbach 
resonance [20]. In order to study this regime quantitatively, it is compulsory to check the 
reliability of the theories adopted in the analysis. Such work has recently been completed 
for the mean-field Gross-Pit aevskii (GP) theories, typically reaching a validity estimate of 
n\a\^ > 10^^, [67,74,75] in combination with an instability criterion for negative scattering 
lengths given by \a\N/bt <~ 0.67, [66]. 

The topic of this thesis is the description of correlations in many-boson systems beyond 
the mean-field. To achieve this, one has to consider not only BEC gas-type states but also 
molecular-type states, since the instability criterion above designates the threshold where the 
latter are formed. In order to simulate a possible Feshbach resonance and break the mean- 
field vahdity region, scattering lengths should be allowed to cover — oo < a < 0. The main 
goal is then to develop the numerical tools needed to understand the nature of interparticle 
relationships in BECs and estimate both the overall importance of such correlation effects 
and the relative importance between the different orders. 

The particular N-body technique chosen for this task is the Stochastic Variational Method 
(SVM). This method provides a solid and arbitrarily improvable variational framework for 
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1. Introduction 



the solution of diverse bound-state problems. A special feature of the SVM is the strategy 
for optimization of a variational trial function by "controlled gambling" . This strategy has 
been proven to be very effective for highly correlated nuclear few-body systems [1]. 

The computational load of the proposed numerical study with the SVM is excessive even 
for the fastest modern computers, and numerical calculations are only feasible for systems 
of three and four bosons. However, since these constitute nontrivial BECs they may work 
as prototype systems in the attempt to describe correlations. In other words, the two key 
questions of the current study are: 

• To what extend does higher-order correlation effects influence the systems of three and 
four trapped bosons ? 

• Can the main conclusions for the four-boson system be generalized to all many-boson 

systems ? 

The remainder of this theoretical thesis seeks the answers to these questions. In chapter 2 
the basic theory needed in a variational treatment of an N-body system is reviewed. The 
variational trial function is a crucial element of this approach. In chapter 3 it is shown 
how to include different levels of correlation explicitly in the functional form of the trial 
funtion. The SVM is introduced in chapter 4, in combination with details of the subsequent 
application to the cases of the He atom and the N-boson systems. Chapter 5 illustrates 
and discusses the numerical results, and the conclusions are collected in chapter 6. The 
derivations of the matrix elements can be found in the appendices. 



1.1 Units and notation 

Where nothing else is indicated, the Atomic Units ^ {rrie = e = ao = h = 1) are used when 
writing results. Moreover, boldface is used for vectors (a) and matrices (A). The length 
of a vector is written \a\ while unit vectors have a hat (a). The elements of vectors and 
matrices are always specified by subscripts (Aij). The elements of a set {A} are sometimes 
most conveniently denoted by superscripts in parenthesis (A^*^^) and sometimes by subscripts 
(Afc). Operators are assigned a wide hat {A). 

With X being an (A^ — 1) x 1 one-column matrix of variables and it's 1 x [N — 1) 
one-row transposed matrix, a quadratic form will be written 

N-l N-1 

tjc .^4cc — ^^^^^ ^^^^^ tJCj^tJC j 

1=1 j=i 

where A is a (N — 1) x {N — 1) symmetric matrix. 

Matrix elements are written in Diracs bra-ket notation 

"'^In this system of units the fundamental electron properties, rest mass m,., elementary charge e, Bohr 
radius ao and angular momentum h are all set equal to one atomic unit (a.u.). This makes the atomic units 
convenient when describing the properties of electrons and atoms or particle systems of similar size. 
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where a denotes all the coordinates of the system and tpi and ipj are square integrable 
functions having a finite scalar product defined by the overlap integral 

Addition of angular momenta is expressed as direct products within square brackets. 
For example the vector- coupling of the angular momenta Ji and J 2, each satisfying the 

eigenrelations Ji<pj,M, = Ji{Ji + ^)(i>J^^■h and Jiz4>j,M, = Miipj^M,, is written 

[0JiAfi ®0J2M2]jM = (^iMi J2M2 I Jl J2 JM) (f)j^Mi J2M2 

Ml M2 

where ( JiMi J2M2I Ji J2 JM) are the Clebsch-Gordan coefficients. 

Computational complexity is discussed in the big oh notation defined [10] 

f{n) = 0[g{n)) f{n) < c * g{n), for all n > rio > and c > 0, 

which means informally that / grows at the same rate as g or slower. 

1.2 Computer programs 

In the course of this work computer programs have been developed in C++ to calculate the 
numerical results (the usage information is listed in appendix lE]): 

• scatlen: Calculates the scattering length for a two-body interaction of identical bosons 

• bee: Calculates the energy for a given state of an N-body system using the Stochastic 
Variational Method. 

The source code for the programs bee and scatlen can be downloaded from my home 
page at: www.phys.au.dk/~hansh. A few examples have been placed in footnotes through- 
out the thesis, indicating the explicit command for the computation of a graph. 
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Chapter 2 

Variational approach to N-body 
problems 

This chapter gives a brief description of how one solves the time independent Schrodinger 
equation for N-body systems using a variational approach. The Hamiltonian is introduced in 
the first section and consists of terms for kinetic energy, two-body interaction and possibly an 
external trapping field. A section then presents the particular two-body interactions applied 
later in the thesis. The most crucial points in the variational method is the introduction 
of a set of relative coordinates and the construction of a flexible trial wave function from 
some appropriate basis of functions. Both points are explained in subsequent sections and 
symmetrization is addressed. Finally, it is shown how the variational theorem defined by the 
trial function reduces to a generalized matrix eigenvalue problem and that accurate results 
can be achieved with basis optimization procedures. 



2.1 Hamiltonian 

In the following, N-body systems of non-relativistic particles are considered, where the ith 
particle has mass m^, charge q, spin Sj, isospin tj and position vector r^. The motion of the 
particles is given by the time-independent Schrodinger equation 

= E^f (2.1.1) 

where the square integrable wave function, \E'(ri, . . . , tat), describes the state of the system 
with the interpretation of as the probabihty density [2]. Most often, this is the common 
starting point of both few-body and many-body treatments. However, the magnitude of 
becomes significant in practice, especially when dealing with identical particles (see below), 
making ab initio restrictions on ^ necessary for many-body systems. 

Assuming the particles move in an external field and that the only particle-particle inter- 
action is through local spin-independent two-body potentials ,Vij, the Hamiltonian becomes 



W , 9 N 

1=1 * t<j 



where = (g^, gf-) is the gradient operator with respect to rj. The subsequent 
separation of the center-of-mass motion from the intrinsic motion allows a translationally 
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invariant description. The following sections introduce the theory necessary to obtain a 
variational solution to ()2.1.H) for the N-body Hamiltonian ()2.1.2|) . 

2.1.1 Identical particles 

Many-body systems often contain a number of identical particles. The indistinguishability 
of identical particles is obviously reflected in the Hamiltonian ()2.1.2j) by the symmetry of the 
operators entering. However, since it is written in first quantization, H does not distinguish 
whether identical particles are bosons or fermions, and therefore this information should be 
added by hand to the wave function, in the form of a definite symmetry. For bosons the 
wave function is required to be even under the interchange of any pair of particle coordinates 
while for fermions it should be odd [2]. Achieving this in many-body problems is only feasible 
with some restrictions on the form of \E' (see section 1^.6. Assuming the proper symmetry 
is given, one may advantageously use a simpler Hamiltonian, given by 

— Hem (2.1.3) 

since all terms in the sums of ()2.1.2j) will contribute the same to the energy. However, 
in the current study, the symmetric "few-body" Hamiltonian in ()2.1.2jl is better suited for 
investigating correlations, and is therefore kept, also for identical particles, in the following. 
The details of applying Hjd with a symmetric (Jastrow-type) wave function for many-boson 
systems can be found in ref. [8], Chap. 2. 



H 



Id 



N -—V, + Ve,t{ri) + -{N-l)V, 



12 



2.2 Two-body interactions 

In sufficiently dilute N-body systems only binary collisions contribute significantly to the to- 
tal energy and three- and more-body interactions can be almost completely ignored (with one 
exception being three-body molecular recombination in atomic gases [71]). The Hamiltonian 
in (j2.1.2j) . introduced as the ab initio starting point of this chapter, takes this simplification 
even further by assuming only spin-independent (central) two-body interactions, Vij. Such 
interactions are sufficient for the calculations of atomic systems and gases of atoms pre- 
sented later. Realistic nuclear models require spin-isospin dependent interactions including 
(at least) three-body terms and are not considered here ^ Now follows a brief introduc- 
tion to the two two-body potentials applied in this thesis and of the concept of the s-wave 
scattering length which is essential in the description of BECs. 



2.2.1 Electrostatic interaction 

The interaction between particles carrying electric charge is the Coulomb force. Like gravita- 
tion, this is a long range one-over-square-distance force, although many orders of magnitude 
stronger. Atomic physics and solid physics, and for that matter the whole of chemistry, can, 
in principle, be determined by this force combined with theories of relativity and quantum 
mechanics [2]. The corresponding interaction potential in SI units is 

V(n2) = (2.2.1) 



^See ref. [8], Chap. 10 for a treatment of the nuclear many-body problem (Argonne/Urbana potentials). 
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Fig. 2.1: Rb-Rb interaction potentials (solid lines) as a function of the atomic separation. Data 
is from Krauss and Stevens [69]. The dashed lines corresponds to the Morse model potential 
[2], V{ri2) = D[(e-°('^i2-ro) _ i)2 _ ^ ^-^j^ ^ 7.87ao, a = OASa^'^ and D = 0.018 a.u., for 
the singlet curve and tq = 11.65ao, a = O.SSoq ^ and D = 0.00093 a.u., for the triplet curve. 
The dash-dotted graph represents the Geltman model for the triplet potential [70], V{ri2) = 
Cg[e-«(ri2-rc)/^6 _ 1/^6^] ^j^j^ ^ 4790 a.u., = 9.7ao and a = 0.90^^ The inset shows the 

details of the triplet potential and the dotted lines indicate the Gaussian models used in this work. 

where qi and q2 are the charges and eo is the permitivity of free space (oq = e^/47reo). The 
non-finite long range character of the Coulomb potential makes solutions of the Schrodinger 
equation difficult in the case of scattering [4] . Moreover, it is not possible to define an s-wave 
scattering length (see below) for this 1/r-asymptotic potential. 

2.2.2 Atom-atom interaction 

The essential property of realistic interatomic interactions is that atoms repel at short dis- 
tances and attract when they are some distance apart. In the following the focus will be on 
the interaction of Rb from the alkali atoms group, since they in particular play a key role in 
experiments on cold atomic gases (and consequently adopted as the default particle in the 
numerical BEC calculations presented later). 

The ground state configuration of alkali atoms has all electrons but one occupying closed 
shells while the remaining valence electron is in an s orbital of a higher shell [6]. In the case 
of ground state collisions, the potential energy depends solely on the internuclear separation 
and the orientation of the two atoms valence electronic spins (sj = 1/2) which couple into 
singlet {S = 0) or triplet {S = 1) configurations, where S = Si-\- S2, [60]. This is illustrated 
for ^^Rb in figure 12.11 The solid curve is based on ab initio calculated data from [69] while 
the dashed and dash-dotted curves correspond to model potentials ^ . For spin-polarized 



^The Morse model potential is an excellent approximation for the short range interaction shape while the 
Geltman model potential (like the Lennard- Jones potential, [2]) has the correct long-range behavior. 
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atoms one may assume that they interact only via the triplet potential shown in details in 
the inset. 

Gaussian model potential 

The key feature of the Rb-Rb potential in the current context, is that it can be assumed 
to have a finite range in terms of scattering. This can be understood from considering the 
quantum mechanical interactions of the Rb2 constituents (a total of 74 identical electrons 
and two nuclei). Clearly, when the nuclei are close enough (about 20A) for the two electron 
clouds to overlap the potential energy will depend greatly on the spin of the valence electrons. 
This is due to the Pauli exclusion principle, since in the triplet configuration the spatial part 
of the electron wave function must be anti-symmetric and so the overlap between electrons 
is minimized in that case. Even at a distance there is residual overlap leading to a long 
range exchange term ^. However, when the nuclei are farther away, the energy due to the 
overlap of electrons decreases exponentially and the interatomic potential is dominated by 
the van der Waals force ^. This dispersion effect can be expanded in a multipole expansion 
such as [68] 

Vd^spirl2) = --^--^--Jo---- (2-2-2) 

'12 '12 '12 

where the leading (van der Waals) term is l/T^g- With such an infinite tail on the potential 
it seems invalid to talk about a finite range and scattering length. Fortunately it can be 
shown that for scattering via power-law potentials, l/r", the decrease is fast enough to be 
considered of finite range provided that n > 3, [6]. 

In this work, the atom-atom interaction is represented by a simple Gaussian model 
potential, V{x) = Vqc"^ , of finite range (the dotted lines in fig. 12.11) . For weak interactions 
in the low-energy (i.e. ultra-cold) limit, the properties of the two-body interaction are 
basically determined by the scattering length, a, introduced in the next section, alone. This 
means, that the exact shape of the potential is insignificant (see e.g. section 13. 3. 3|) and that 
the apparent lack of a hard core at small ri2 in the Gaussian model is acceptable. More 
details of the Rb-Rb interaction can be found in [60]. 

2.2.3 Scattering length for finite interactions 

The following is a very brief account of basic scattering theory which can be found in 
details in Refs. [2,3]. Consider the situation of two isolated particles with masses mi and 
m2 that interact via a central potential ^^(^12), where ri2 = |ri — is the interparticle 
distance. Further assume that the interaction vanishes rapidly (faster than oc rfg^) for large 
separations, i.e. V{ri2) —>■ for ri2 —>■ 00. As outlined in section 12.31 the motion of the 
particles separates into the trivial center-of-mass motion and the relative motion described 
by a single coordinate wave function, \E'(a;), satisfying the Schrodinger equation 



-2 



^(a;) = E^(a;) (2.2.3) 



•^This term has the form, Vex{ri2) — ±Ar^2e~'^^'^^^ , also found in the Morse model potential. 
^As the electrons move, small fluctuations occur in the charge density surrounding each atom so, in turn, 
one atom can polarize the other setting up an momentary dipole moment which then attracts the first. 



2. 3 Relative coordinates 
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where x = Vu and /i = 17111712/ {mi + 1712) is the reduced mass. Solutions with E < 
correspond to bound states of the potential. Scattering is described by the Lippmann- 
Schwinger solutions [3] with positive energy E = h'^k'^/2fi, and the asymptotic form 

^ [e-- + m, kf-^] (2.2 A) 

corresponding to the sum of an incoming plane wave with relative momentum hk and a 
scattered spherical wave (i.e. the (+) superscript) with amplitude 

f{k',k) ^ — ^ y dx'j^^V{x')^'^:\x') (2.2.5) 

For a spherically symmetric potential the scattering amplitude only depends on the angle, 6, 
between the relative momentum of the particles before and after the scattering, f{k' , k) = 
f{k, 6). In the low energy limit, 0, where isotropic s-wave scattering is dominant ^, the 
scattering amplitude approaches a constant, /(O, 0) = —a, and the wave function reduces to 

^W(a;)'=°^W(a;) = l-^ (2.2.6) 

The constant a is the s-wave scattering length and can thus be determined as the interception 
of the asymptotic wave function and the x axis, that is '^{a) =0 for the zero energy solution 
to the Schrodinger equation. Figure 17!^ ^ demonstrates this in the case where two identical 
particles are interacting via a finite Gaussian potential, V{x) = VqC"^ . The scattering 
length is always positive and finite for repulsive interactions, Vq > 0, while for attractive 
interactions, Vq < 0, it can be both negative and positive and becomes divergent when 
changing sign. This behavior (zero-energy resonance) occurs each time the potential is just 
deep enough to support a new bound state. 



2.3 Relative coordinates 

The most convenient way to remove the center-of-mass motion is to express the Hamiltonian 
in terms of relative coordinates that do not change when the system moves or rotates as a 
whole. The obvious choice is the scalar interparticle distances 

r^, = \r^-rJ\, z^j = l,...,iV (2.3.1) 

where is the position of the ith particle, giving A^(A^ — l)/2 relative but dependent 
coordinates. Choosing one relative coordinate, ri2, in the two-body case is trivial. In three- 
body problems the truly independent positive perimetric coordinates 

Ui = ^{rik + rij -rjk), i j k = (1,2,3), (2.3.2) 

partial wave expansion of '^{x) in Legendre polynomials, Pi{cos9), gives a radial Schrodinger equation 
where the effective potential includes a centrifugal barrier, i.e. the term h'^l{l + l)/(2/ia;^), [4]. Thus waves 
with energies much lower than this barrier are simply reflected leaving only the s-wave (/ = 0) contribution. 

^The numerical data for fig. 12. 21 is produced by calling the program scatlen. Graph a) is generated with 
the command lines: scatlen -VO 2.11e-7 -compare, scatlen -VO 2.11e-7 -printwave and scatlen 
-VO 2. lle-7 -printpot. For b), c) and d) just change 2. lle-7 to the corresponding values of Vq. 
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Fig. 2.2: Scattering lengths, a, and wave functions, ^(2;), for Gaussian potentials Vix) = Voe~^ ' , 
with b = 18.9 fixed and different strengths: a) repulsive, Vq = 2.11 x 10"^; b) weak attractive, 
Vq = —2.11 X 10~^; c) more attractive at bound state threshold, Vq = —4.743 x 10~^; d) strong 
attractive with 4 bound states, Vq = —2.11 x 10~^. The functions have been scaled to fit [—1; 1]. 
Other numbers are in atomic units. 



are often preferred as this simplifies integral evaluations over the coordinates (used in ap- 
pendix [HiBI). Since there is only 3A^ — 6 internal space degrees-of-freedom in an A^-body 
problem (for N > 2), the set of scalar relative coordinates include unnecessary extra coor- 
dinates when N{N — l)/2 > 3A^ — 6 -v^ > 4. This complicates the use of interparticle 
coordinates in many-body systems significantly [46]. 

A different approach, also convenient if A^ > 4, is to introduce a set of relative vector 
coordinates = {xi,X2, . . . ,Xi^_i) and the explicit center-of-mass coordinate x^. They 
are related to the single-particle coordinates = (ri, r2, r^r) by a linear transformation 



N 

Xi = Y,Uijrj, t = l,...,N (2.3.3) 

i=i 

written in matrix form as a; = Ur, where t/ is a suitable N x N transformation matrix . 
A widely used choice for x, which is also employed for the many-body problems considered 
in this thesis, is the Jacobi coordinate set [8] defined by 

Xi = ^i{Ci~ri+i), i = l,...,N - 1 (2 3 4:) 



where Ci is the center-of-mass and Uj = the reduced mass of the first i particles. 

' mi2...i+i ^ 

^Obviously one can readily generalize the scalar definition H2.3.1|) to an equivalent vector description with 
= ('"12, '''13, I'lN, ''"23, •■, '"2Af J f(N-i)Ni Xm), given by an appropriate linear M x N transformation of 
the single-particle coordinates. 



2. 3 Relative coordinates 
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The corresponding transformation matrix is 



/i2 — 

mi2 



C/.i 



mi 



V 



mi2...]v-i 
mi 

mi2...]v 



V/^Af-l 



m2 



mi2...]v-i 
m2 
mi2...]v 



ms 



mi2...jv_i 
ms 

"1.12. ..AT 








\ 



mi2...jv / 



where the short notation means mu-.i = mi + m2 + ■ ■ -mj, making mi2---N the total mass 
of the system. A specific set of Jacobi coordinates, {a:;,}, are related to the interparticle 
distances, rij, and the hyperradius, p, through the relation [64] 

N N N-1 

i<j i=l i=l 

Moreover, since the first Jacobi coordinate, Xi, of the set defined in ()2.3.5p is equal to 
.y//Ii(ri — = y/JIiri2, the two-body interaction, Vu, depending on this interparticle- 
distance, has a simple form V{xi/ ^/JIl). By permutations of the particle labels (1, 2, ■ ■ ■ ,N), 
and the corresponding columns in ()2.3.5|) . one can generate different Jacobi coordinate sets 
x^P\ known as arrangements or partitions, where the first Jacobi coordinate is y/JJiVij and 
the form of Vij is simple. This is an advantage when trying to obtain analytical expressions 
for the integrals (see ()2.4.6p and ()2.4.7p derived in the next section) that have to be evaluated 
in the variational method. 



2.3.1 Separation of the center-of-mass 

The many-body Hamiltonian ()2.1.2|) written in terms of relative coordinates separates into 
a translationally invariant part and a part involving only the center-of-mass coordinate. 
Corresponding to the change of coordinates ()2.3.3|) the single-particle gradient operators 

V =(Vi,V2,..,V7v), entering the kinetic energy part, are transformed by 

V = C/^V^ (2.3.7) 
and since the transformation matrix U can be assumed to satisfy the relations U at,- = — — — 

mi2...]v 

and Uij = 6Ni [1], easily verified for Uj, one has 



1=1 * i=i * fc=i 1=1 

^ ,^ rr ^ zmi2...N 

k=l 1=1 1=1 
^2 N-1 N-1 

= -y5Z5Z^'^''^-^'^-'+^-™ (2.3.8) 

k=l 1=1 
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where Tcm = — = —75—^- — is the center-of-mass kinetic energy ^, and 

2mi2...jv 2mi2...jv ^JV ' 

^ TI TJ 

A« = V^^, fc,/ = l,...,iV-l (2.3.9) 

1=1 

The external field potential, Vexti separates in a similar way when applying transformation 
()2.3.3|) . see section 031 Thus the Hamiltonian ()2.1.2|) can also be expressed in terms of the 
quadratic form 

N-l N-l 
k=l 1=1 

as 

H = -y V^AV, + ye.t{Xi) + ^^J- + Hem (2.3.10) 

1=1 i<j 

with Vext and Vij depending on the relative coordinates {xi, X2, ■ ■ . , x^-i} and Hem = 
Tcm + Vext{R)- Explicit insertion of the Jacobi transformation matrix, Uj, in the expression 
fl2.3.9p for A, produces the important result A = /. This means, that in the specific case 

of the Jacobi coordinates the quadratic form is dissolved and only terms of the Laplacians, 

^ 2 

V , have to be considered in a calculation of the kinetic energy. 



2.4 Matrix representation 

The general stationary-state solution, \E', to the Schrodinger equation p.l.lj) will be a linear 
superposition of the eigenfunctions \E'n of H: 

00 

* = ^a„^„, (2.4.1) 

n=l 

where satisfies 

^*„ = E„^„, n = l,2,.... (2.4.2) 

If the Hamiltonian (j2.3.1(Jj) is Hermitian, the eigenvalues En are real and the eigenfunctions, 
"^n, called the energy eigenstates, form a complete and orthogonal set {^E^n}, [2]- With focus 
on bound-state solutions, in particular the ground state and lowest excited states, one has 
to find the lowest discrete energies. En, and corresponding eigenstates, "^n, from eq. ()2.4.2|) . 
Unfortunately, except for the two-body cases, the explicit form of the eigenstates is not 
known a priori, making it impossible to solve the eigenvalue problem analytically. 

The best alternative is to consider a finite set of known functions {ipi, -^2, • • • , 4'k}, that 
are linearly independent and possibly non-orthogonal. A general function in the space Vic 
spanned by this set can be written 

K 

^ = 5^c,V^^ (2.4.3) 

i=l 

^Corresponding to the total momentum p^^i — ~ = ~*^X^jLi ^Nj^x, — —i-fi^xN- 



2.5 The linear variational method 
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The state vector c uniquely defines ^ in the function space Vic- Inserting this form into the 
Schrodinger equation gives 

K 

Y,{H - e)c,ijj = 0, (2.4.4) 
i=i 

which is the eigenvalue problem for H inside Vfc- From this restricted problem one can 
determine eigenvalues €1,62, ... ,eK and corresponding state vectors c^^\ c^'^\ . . . , c^^) that 
are approximations to the exact solution ()2.4.2|) . As will be discussed in the next section it 
is actually the best solution within Vn from a variational standpoint. 

Equation ()2.4.4j) is conveniently expressed in a matrix representation by multiplying from 
the left by ip* and integrating over all coordinates ol that ipi depend on, giving 

K 

Y^iHi, - eS,,)c, =0, t = l,2,...,K. (2.4.5) 

where 

= {iPi\H\^j) = J ri{<^)Hij,{cc)dcc (2.4.6) 
are the elements of the K x K Hamiltonian matrix and 

S,, = (V-^,) = / V>*(a)V',(a)rfa (2.4.7) 



are the elements of the K x K overlap matrix. In this way, the problem of determining the 
eigenvalues of the operator H in ()2.4.4|) has been transformed to a generalized eigenvalue 
problem of square matrices ^ , most elegantly written 

He = eSc. (2.4.8) 

If the basis functions are chosen to be orthogonal, S becomes the identity matrix, and 
equation (j2.4.8p reduces to a standard eigenvalue problem. 



2.5 The linear variational method 

In this thesis, a variational method is used to obtain the approximate bound-state energies 
and wave functions of systems described by the N-body Hamiltonian (|2.1.2j) . It is linked 
to the fact that an arbitrary wave function corresponds to an energy higher or equal to the 
true ground state energy. The so-called variational theorem states 

for any square-integrable ^. The functional £ is the expectation value of H and the equality 
holds only if \1/ is the ground state of H with the eigenvalue Ei. A proof of this theorem is 
elementary and can be found in textbooks on quantum mechanics, e.g. [2] p. 116. 

^The connection between the algebra of linear operators and square matrices is quite fundamental, see [7] 
sec. 5.10. 
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The variational theorem is the basis of the widely used Rayleigh-Ritz variational method. 
The idea behind this method is to choose a trial function as \1/ that depends on a number of 
variational parameters. Evaluating the expectation value S yields a function of these param- 
eters, and by minimizing with respect to the parameters, one obtains the best approximation 
to El that the explicit form of \1/ allows. 

The linear variational method is a variant of the Rayleigh-Ritz method in which one 
works with trial functions of the form ()2.4.3|) . described in the previous section. The expan- 
sion coefficients = {ci, C2, . . . , ck} serve as the variational parameters and minimization 
consists of demanding that 

1 = and ^ = 0, (2.5.2) 

for all i = 1,2, . . . , K. Considering first the latter condition in ()2.5.2|) by differentiating S 
with respect to c* gives 



dc* (^1^)2 



For this expression to vanish the numerator must be zero. Introducing the expansion ()2.4.3p 
as \1/ in the matrix element expressions, one has 

K K K K 

mH\^) = 5^5^c*c,(^.|^|^,) = 5^5^c*c,i/., (2.5.3) 

i=l j=l i=l j=l 

and 

K K K K 

= 5^5^c:c,(^#,) = 5^5^c*c,^., (2.5.4) 

i=l j=l i=l j=l 

where Hij and Sij are again the elements of the Hamiltonian and overlap matrices, and the 
condition ^ = reduces to 

d d 

g^imm - ^^(^1^) = - ^B^^y^ = 0' ^ = 1, 2, . . . , 



The condition ^ = will reduce to the complex conjugate of this equation [7] and hence 
gives no new information. Thus imposing the conditions ()2.5.2|) on the expectation value £ 
has produced precisely the same matrix eigenvalue problem 

He = eSc (2.5.5) 

that was derived in the previous section. This is an important connection and means that 
by choosing an arbitrary basis {iIji,iIj2, ■ ■ ■ , "^Pk} of known functions and solving eq. ()2.5.5|1 . 
also called the secular equation, one will in fact get the best approximate solution inside 
Vk- One could say, that the minimization with respect to the parameters {ci, C2, . . . , ck} is 
implicit in the solution. Of course, the functions ipi can be made dependent on additional 
nonlinear variational parameters, giving further fiexibility to the trial function. 



2. 6 Basis functions 
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The variational theorem implies that the lowest of the eigenvalues determined by eq. 
fj2.5.5|) will be an upper bound to the real ground state energy Ei. In turn, all the eigenvalues 
ei,i = 1,2, . . . , K, are upper bounds to eigenstate energies of the full Hamiltonian (see [1], 
theorem 3.3). Arranging in increasing order the K eigenvalues ei < £2 < ■ ■ ■ < of the 
truncated problem and the discrete eigenvalues Ei < E2 < . . . of the full problem (j2.4.2p . 
it can be shown that 



Expanding Vk by increasing the number of functions in the basis will bring fl2.5.5p closer to 
the full Hilbert space problem, and obviously improve on the approximate eigenvalues by 
lowering them towards the exact values Ei 

2.6 Basis functions 

A crucial point when using the linear variational method is the choice of basis functions. An 
expansion in the basis should give a good representation of the physical shape of the wave 
function for the quantum system in question. It is important, in general, that two basic 
requirements are satisfied: 

• The basis should form a complete set so that the result obtained by a systematic 
increase of the number of basic functions will converge to the exact eigenvalue. 

• Furthermore, all matrix elements should be analytically calculable, for the variational 
approach to be practical. 

In addition, to solve an N-particle problem accurately and with a high convergence rate, 
the trial function, \E' = J2f=i^i'^iy hence the basis functions, {^ipi}, should describe 
the correlation between the particles well, have the proper symmetry and encompass the 
appropriate degrees of freedom, e.g. orbital and spin angular momenta. In this thesis, it is 
assumed that a basis function with total angular momentum J and projection M, can be 
written in the form [1,28] 



where P is a sum of permutation operators ensuring the proper symmetry, describes 
the spatial dependence, 9lml{^) specifies the orbital motion with definite angular momen- 
tum L, and Xsms Vtmt the spin and isospin parts Since all correlations between 
particles in the systems treated later are state- independent, they can be fully represented in 
the spatial part, ^^.{x), of ipk, as discussed in detail in chapter|21 More elaborate descriptions 
of the angular momentum parts can be found in appendix 1X1 

2.6.1 Symmetry 

From experiments it is known that particles of zero or integral spin, such as the photon and 
^He, are bosons. Particles with half-integral spin values, such as electrons and nucleons, 

^°This can be explicitly proven, see [1], p. 27. 

^^Parts for other degrees of freedom, like color and flavor, can be added correspondingly. 



El < ei, E2 < 62, ■ . ■ , Ek < eK 
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are fermions. Atoms constitute bosons when they contain an equal number of nucleons, 
otherwise fermions. When a system consists of a number of identical and indistinguishable 
particles the wave function must have the proper symmetry with respect to any interchange 
of the space and spin coordinates of the identical particles. The proper symmetry of the trial 
wave function can be achieved by operating on the basis functions with the operator [8]; 



using ap = 1 for identical bosons and ap = (—1)^, where p=0,l is the parity of the permu- 
tation, P, for identical fermions. Here, the permutation operator, P, permutes the variable 
indices (1,2,..., N) of identical particles to (pi,P2, • • • ,Pn) and the summation over P in- 
cludes all necessary permutations. Thus V corresponds to a symmetrizer {S) for bosons and 
an antisymmetrizer (A) for fermions. 

Obviously, the permutation operator, V, commutes with the symmetric many-body H 
in ()2.1.2j) . but, except for the case N = 2, the A^! different permutation operators do not 
commute among themselves. This means, that an eigenfunction of H is not necessarily an 
eigenfunction of all P. Only a totally symmetric eigenfunction, \E'5, or a totally antisymmetric 
eigenfunction, can be common eigenfunctions of H and all P. The two types of wave 
functions, "^s and "^a, are thought to be sufficient to describe all systems of identical particles 

Accordingly, if a system is composed of different kinds of identical particles, its wave 
function must be separately totally symmetric (bosons) or totally antisymmetric (fermions) 
with respect to permutations of each kind of identical particles, [2]. 

2.7 Basis optimization 

The most direct approach to a variational solution of a quantum mechanical bound-state 
problem is to solve the secular equation defined by a basis of functions, {ipi,4'2, ■ ■ ■ yi'K}, 
containing no nonlinear parameters. Two main steps are involved: the calculation of the 
overlap and Hamiltonian matrix elements and the solution of the generalized eigenvalue 
problem. To achieve the desired accuracy one only needs to add many (linearly independent) 
functions to the basis. Unfortunately, calculating all matrix elements takes time (9(A'^) on 
a computer and solving eigenvalue equations is an 0(i^'^)-procedure, [39]. Any variational 
approach is thus only feasible on a basis of reasonable size consisting of functions that 
allow fast matrix element evaluation. The direct method in particular suffers from this 
limitation since the convergence is often slow increasing the demand for a large number of 
basic functions, [43]. 

Another approach, designed to avoid a huge basis dimension, is basis optimization. The 
idea is to only select the specific basis functions that give good results. To this end, the shape 
of the basis functions is made dependent on nonlinear parameters, which, in effect, determine 
how well the variational function space, Vk, contain the true eigenfunction. An optimal basis 
of definite size, K, can be established by minimizing the variational energy function with 
respect to these parameters. However, although numerous elaborate methods are available 
for multidimensional function minimization (see [40]), the optimization of the nonlinear 

^^This is the so-called symmetrization postulate, [2]. 




(2.6.2) 



2. 7 Basis optimization 
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Fig. 2.3: The control flow diagrams for two common basis optimization strategies: (left) Optimizing 
while increasing one basis at a time and (right) optimizing the parameter generation or intervals. 



parameters in a trial wave function is by no means a trivial task. In fact, computational 
complexity studies show that the general problem takes time exponential in the number 
of parameters and it is therefore rated as intractable, i.e. not amenable to a practically 
efficient solution. This means that the sheer number of nonlinear parameters quickly becomes 
the bottle-neck in basis optimization. 

Different strategies have been employed for basis optimization in few-body problems 
related to the following scenarios: 

• In few-body problems the number of nonlinear parameters needed in each basis function 
is reasonably low, and one might be able to perform a full simultaneous optimization 
of the basis parameters. The strategies employed can be placed into two categories: 
the deterministic and the stochastic. The former is based on stepping or gradient 
strategies (e.g. the conjugate-gradient method, see [33]) and are sensitive to a given 
starting point, always reaching the same minimum from the same initial condition. The 
solution in these cases may not be the global minimum sought but a local minimum. 
Convergence depend heavily on the initial guesses for the parameters. A stochastic 
minimization tends to convergence much slower but eliminates the risk of ending up 
in a local minimum [1]. Combining analytical gradient and stochastic techniques in a 
mixed approach, as proposed recently in [34], seems to be very economical. 

• When the above approach is too time-consuming, which is most often the case, one 
can use grid methods to reduce the number of parameters to a smaller number of 

^^Vavasis [11] reports the worst-case complexity of minimizing a Lipschitz constant function, 
f{xi,X2, --jXa), in a box to be 0{{-^y'), where L is the Lipschitz constant. 
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tempering parameters, either by fixing parameters through a geometric progression [42] 
or pseudo-randomly [37]. In the latter case this corresponds to optimizing the hmits of 
the intervals from which the pseudo-random numbers are selected [45] . The drawback 
is, that parameters may be assigned values disregarding whether the corresponding 
basis functions contribute to the solution or not. 

• Alternatively, a partial optimization can be performed where only a few parameters are 
optimized at a time and all others fixed. In particular, if only one specific basis function 
is optimized, then only one row of the (symmetric) Hamiltonian and overlap matrices 
are affected. Even the consecutive solving of the full generalized eigenvalue problem 
can be avoided in the optimization procedure. This is employed in the Stochastic 
Variational Method and, as described in the chapter HI takes only a fraction of the 
computational time of a full optimization. 

Fig. 12.31 shows two control flow structures that can be used with the optimization strategies 
discussed here. The left diagram corresponds to the case where one basis function is opti- 
mized and added to the basis at a time and the right diagram is for a procedure where the 
entire trial function is constructed and subsequently optimized. The dashed boxes desig- 
nate a (possibly complicated) optimization method in which place the SVM trial and error 
technique will be considered in chapter HJ 



Chapter 3 

Correlations in many-body systems 



The basic theory necessary for treating particle correlations in N-body systems with a varia- 
tional approach is presented in this chapter. The main goal is to develop several descriptions, 
each representing a different level of correlation, and allow for a direct comparison of the 
corresponding correlation energies. This is based on the vital assumption, that correlations 
can be explicitly included in a description, by embedding them, ab initio, in the form of the 
variational trial function. To this end, respective sections treat first the uncorrclatcd Hartree 
trial function used within the Hartree- Fock theory, then the effective interaction (mean-field) 
approach based on the pseudopotential approximation and last an explicitly correlated trial 
function designed to handle two-body, three-body and higher-order correlations. To begin 
with, however, a short remark about correlation as a concept. 

3.1 Defining correlations 

Since there are various definitions of the term correlation available in the physics literature, 
it is appropriate to define the concept clearly before deriving a theoretical description. In the 
dictionary, correlation is explained as "a shared relationship" or "causal connection" . Within 
the physics context of N-body systems, correlation correspondingly designates the possibly 
complex intcrparticle relationship among the particles. However, in some textbooks, the 
energy connected with such correlated behavior, Ecorr, is defined as the difference between 
the energy of the similar non-interacting system and the exact measured or calculated energy 
[4]. In other theoretic areas, hke the atomic Hartree- Fock theory, the correlation energy is 
regarded as the difference between the energy obtained with an independent particle model 
based on the Hartree product wave function (see below) and the exact energy [2,7]. Although 
both interpretations have valid argumentation, they are also very distinct on the important 
question of what defines an uncorrelated system. 

Here, and in the remainder of this thesis, the following definition is adopted: 

In an N-body system, where the interaction between the particles is state-independent, 
the (inherently) correlated motion of the particles can be represented by a wave func- 
tion of the form [8] 

^'(ri, r2, . . . Tn) = F(ri, rs, . . . r7v)$(r-i, ra, . . . Vm) (3.1.1) 

where F is a correlation factor and $ is an uncorrelated wave function corresponding 
to a system of independent particles. The specific correlation energy included in such 
a representation, Ecorr, is defined by 
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N 

Ecorr = \Einteraction\ = I < ^| Vjj \ ^ > | (3.1.2) 

i<j 

where Vij is the two-body interaction potential. 

With this definition, the energy reference point corresponding to an uncorrelated system 
is given by the energy of the non-interacting system. Contrary to other interpretations, 
this allows even a mean-field theory with F = 1 to represent correlation energy, since any 
type of interparticle interaction is tantamount to correlation effects. The belief of the writer 
is, that such a standpoint is more true to the "civil" perception of the word correlation, 
and in any case, advantageous in the current context, because the primary aim here is to 
compare different levels of correlation where the widely used mean-field approach is just one 
candidate. 



3.2 Hartree-Fock mean-field description 

The independent particle model, originally formulated by Hartree in 1928 [58] and generalized 
with symmetry by Fock and Slater [59], is based on the ansatz that a many-body wave 
function can be written as a properly symmetrized product of orthogonal single-particle 
states, given by 

N 

'^HFiri,r2,...rN) = V'^Hiri,r2,...rN) = VY[Hr^) (3.2.1) 

i=l 

where V is defined in ()2.6.2|) as the symmetrizer for bosons and the antisymmetrizer for 
fermions. In the Hartree-Fock method [59] this form of wave function is applied with the 
variational theorem in (j2.5.H) . by demanding that the variation of the energy functional value 
is zero, i.e 6£hf[^ hp] = S < hfIH]"^ hf >= 0. To derive the theory of this method in 
the case where H is the non-relativistic N-body Hamiltonian (|2.1.2|) . one may conveniently 
rewrite H as 

N N 1 2 

H = J2^^ + J2^^^^ where /i, = --— v' + Kxt(r'.) (3.2.2) 

1=1 i<j 

having the first term explicitly given by a sum of N identical one-body Hamiltonians, h^. 
Taking into account that H is invariant under the permutation of particle coordinates, i.e. 
[H,V] = 0, and =V hj definition, the expectation value of hi is simply 

< ^HFlhil'^HF >=< '^H\hiV\'^H >=< (f)i\hi\^i >= J dr(j)*{r)hi(j)i{r) (3.2.3) 

assuming the single-particle states are orthonormal, < >= 6ij. Using the same argu- 

ments the expectation value of the two-body interaction, Vij, becomes 

< '^HF\Vij\'^HF > = < "i/HlVijVl^H > = < MjlVijlMj ± > (3.2.4) 



3.2 Hartree-Fock mean- field description 
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where ± indicates a + for bosons and a — for fermions. Here a two-particle matrix element 
involves a double integral over the coordinates of both particles 



< (f),(j),\V,,\M, >= / / drdr'(j)*ir)(f)*{r')V{r - r')(f),{r)(f),{r') (3.2.5) 




Proceeding by taking the variation of Ehf with respect to the single-particle states, (pi, while 
imposing the orthonormality constraints on the 0j's by introducing (diagonal) Lagrange 
multipliers, £i, yields 

N 

(5EHF-5^^^5<0^|0i>=O (3.2.6) 

i=l 

After some algebra (see ref. [2]) this variation leads to the N Hartree-Fock integro- differential 
equations for the single-particle wave functions 

hMr) + Vgp{r)<P,{r) + j dr'Vi^{r, r')<j),{r') = £,<j).,{r) (3.2.7) 

where the direct potential is 

N 



VSpir) = J2f dr'<P*{r')V{r - r')<p,{r') (3.2.8) 

j=i 



and the exchange potential is 

N 

Vf|;(r,rO = ± 0*(r')V(r - r')0,(r) (3.2.9) 
i=i 

with + for bosons and — for fermions. Adding the requirement of self-consistency between 
the approximate individual single-particle states, 4>i{r), and the variational interaction po- 
tential, V{r — r'), the equations ()3.2.7|) - ()3.2.9|) can be solved with a simple iterative proce- 
dure [2]. 



3.2.1 Correlations in the Hartree-Fock description 

The derived Hartree-Fock equations provide useful physical insight. First of all, they show 
that if the N-body wave function is approximated by the single-particle product ()3.2.1|) . 
the corresponding variational solution describes a model where each particle moves in an 
effective potential generated by the other — 1 particles (i.e. a mean-field). A further 
striking feature of the integro-differential equations is that they involve the joint probability 
for finding particles in states i and j at points r and r'. This obviously imposes a relationship 
among the coordinates of the particles which indicates they are to some degree correlated ^. 

The direct term represents the average potential due to the local presence of the other 
particles. The exchange term takes into account the symmetry effects from exchanging 
particles and indicates that the effective single-particle potential is both state dependent 
and nonlocal. Determining one 0i(r) requires the states for all other particles throughout 

^One may note that the correlations induced by the exchange term are repulsive for fermions (on a range 
comparable to the size of the system) and corresponds to the Pauli blocking effect. 
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the system as well as all other r'. This means that the independent particle approximation 
does in fact not entirely neglect particle-particle correlations. Rather it assumes that most 
of their important effects can be taken into account with a sufficiently clever (variational) 
choice of the two-body interaction potential form Vij. As explained in the next section the 
optimal potential is not the exact particle-particle interaction. 

It is clear, that in the first quantized Hartree-Fock derivation given above, the only dis- 
tinction made between bosons and fermions is the definite symmetry of the wave function. 
This seems only to have minor implications given by a sign in the exchange potential. How- 
ever, at the ultra-low temperature quantum level this difference in the exchange correlations 
of bosons and fermions becomes very pronounced. While bosons eagerly fall into a single 
quantum state to form a Bose-Einstein condensate fermions tend to fill energy states from 
the lowest up, with one particle per quantum state. To exemplify and complete the Hartree- 
Fock description, an expression for the ground state energy of the N-boson system is now 
derived, since this is the case of interest later. The corresponding derivation for systems of 
identical fermions, which are not considered further here, is briefly addressed in appendix iBl 



3.2.2 Ground state of identical bosons 

Bosons in a many-body system obey Bose-statistics with no restrictions on the allowed 
quantum states. The ground state for identical bosons will then have all particles occupying 
the lowest orbital, (j)i{ri) = ipo{ri), and the symmetric Hartree wave function 

^Sf = ^0(^1)^0(^2) . . . MrN) (3.2.10) 

is appropriate as the starting point of the Hartree-Fock method. The spin part of the wave 
function, left out here, is similarly a product of single-boson spin functions but otherwise does 
not enter the calculation. Since the permutation operator, V, is superfluous in the ground 
state derivation, the exchange term in the integro-differential equations (j3.2.7|) vanishes. 
Removing the self- interaction contribution {i = j) from the direct term the ground state 
Hartree-Fock equations for identical bosons (mj = m) is then 

- I^V' + Kxt(r) + (iV - 1) y dr'roir')Vir - r')Mr')]Mr) = l^Mr) (3.2.11) 

where m is the mass and n corresponds to the chemical potential encountered in the Bo- 
goliubov theory. The total ground state energy of the system of N identical bosons becomes 

= Nfi- < ^oMVijl^Po^o > (3.2.12) 

where the second term is due to double counting (see app. El- Assuming that the system is 
sufficiently cold and dense for the single-boson wave functions, ipo{f), to overlap, the Hartree 
wave function, \E'^'', corresponds to a condensed state, as explained previously. Then E^^^p 
is the Hartree-Fock approximation to the BEC energy and apparently takes the boson- 
boson interactions into account. However, another consequence of interactions is collisional 
excitations, where the bosons are scattered in and out of their single-particle states, leading 
to quantum depletion of the lowest state even at T = 0. This is entirely neglected in the 
derivation of (j!mT|l and (jlTTT^ . 



3.3 Pseudopotential mean-field description 
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3.3 Pseudopotential mean-field description 



A wave function in the Hartree form ()3.2.H1 is explicitly uncorrelated. Consequently, it is 
clear that the Hartree-Fock description in the previous section ignores real dynamical corre- 
lations although incorporating exchange and mean-field effects. In the dilute limit, however, 
it is possible (and in some cases imperative, see below) to introduce some extend of (short 
range) correlation effects in the variational solution by adopting an effective interaction po- 
tential instead of the exact Vij. The details of this important refinement of the Hartree-Fock 
theory are presented in the following. 

As a first approximation one must assume that the system consists of weakly interacting 
particles (what exactly constitutes a weak interaction will be addressed in section I3.3.4|) . 
When this is the case, the interaction potential is so short that the single-particle wave 
functions do not vary over the interaction region ^. Then one can rewrite, as 



dr'(l)*{r')V{r - r')0j(r') ^ \(pjir)\^ j dr'Vir - r') 
simplifying the intergro-differential Hartree-Fock equations ()3.2.7|) to 



hi + VMpir 



Sidi(r) 



(3.3.1) 



(3.3.2) 



where the mean-field potential is 



N 

VMF{r) = J2(^±6,,M,{r)\' / dr'Vir -r') 



(3.3.3) 



where -|- is for bosons and — is for fermions. Thus, in the bosonic ground state described 
above, the exchange term simply doubles the direct term. For identical fermions in equivalent 
single-particle states the terms cancel instead as expected due to Pauli exclusion. 



3.3.1 Effective interactions 

Intuitively the mean- field interaction, Vmf, should be mediated by the elastic collisions in 
the system. As mentioned, only two-body scattering described by the two-body potential, 
Vij, are significant in a dilute system. However, there are several reasons not to use the exact 
two-body interaction potential in the Hartree-Fock approach. First of all, it is quite difficult 
to determine the exact potential precisely, and a small error in the shape of Vij may produce 
a large error in other results, e.g. in the scattering length, a. Secondly, the exact potential 
is very deep and supports many bound states. Such strong interactions cannot be treated 
within the weak field assumption (see below). Finally, the hard-core repulsion at short 
distances makes the evaluation of the mean-field integral (j3.3.3j) somewhat troublesome. For 
an extreme example, consider the hard-sphere potential {V is infinite if (r — r') < Tc, zero if 
(r — r') > Tc). Obviously this potential causes Va//f(^*) to be infinite for any nonzero value 
of (piir), regardless of the size of the hard-sphere radius Vc- This is unreasonable since even 
in the limit of an infinitesimal radius, the contribution would still be infinite. 

^In the independent particle model weakly interacting particles are roughly free particles given by plane 
waves, (/>i(r) = (27r)~^/^e*'" ''. Thus the approximation (f>i{r) « (t>i{r') amounts to the requirement that the 
thermal de Broglie wavelength, At = (^Tsr)^^^; is much larger than the range of V{i — r'). 
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The explanation for the discrepancy of the hard-sphere example is that in the independent 
particle approximation, no dynamical correlations between individual particles are allowed. 
In reality, there would be no particles closer to each other than the radius Vc in the hard- 
sphere scattering example. However, there is no way for the naive Hartree-Fock theory to 
account for this. Simply neglecting the hard-sphere interaction, as there are no particles 
that close anyway, is not sensible either since this would allow the single-particle wave 
functions of two neighboring particles to overlap inside the hard-core radius ^. The solution 
is instead to replace the exact interaction potential by a model potential that (1) has the 
same scattering properties at low energies, i.e. is the same scattering length and (2) will 
work in the independent particle approximation. To some extend, the short wave length 
components of the wave function that reflect the dynamic correlations between particles are 
then implicitly taken into account. This implies that the Born approximation in the case 
of scattering (see below) and the Hartree-Fock method for calculating bound state energies 
give better results provided that the simple effective interaction is used rather than the real 
one. To exactly what extend such a mean-field approach succeeds in including correlations 
is somewhat clarified in chapter [SI 

3.3.2 The pseudopotential approximation 

The model potential satisfying the two requirements stated above with the minimal number 
of parameters (one!) is the zero-range pseudopotential initially introduced by Fermi [72] and 



where the coupling constant, g = ^^a, is directly proportional to the s-wave scattering 
length, a. It is valid for dilute systems (typically stated as n\a\^ -C 1, where n is the 
characteristic density) at low energies, although making g energy or density dependent can 
extend the validity region [64,75]. The pseudopotential involves a Dirac ^-function and a 
regularizing operator, g^ri2, that removes a possible divergence of the wave function at 
ri2 = 0. When the wave function is regular at ru = the regularizing operator has no effect 
and the pseudopotential can be viewed as a mere contact potential, V^r^) = g5{r 12). The 
widely used Gross-Pitaevskii equation ^ corresponds to a mean-field approach for systems of 
bosons based on the pseudopotential effective interaction. 

In the low energy limit the pseudopotential reproduces the scattered wave function of the 
exact two-body potential asymptotically and gives the correct scattering length. However, 
the possible change of the wave function inside the (finite) interaction range is effectively 
ignored. This is also known as the shape-independent approximation, [63]. As demonstrated 
above with the hard-sphere potential example, it is the large repulsive core of the exact inter- 

■^The single-particle wave functions need to bend away in the forbidden hard-core region so that the 
resulting curvature of such </)i's contribute correctly to the kinetic energy of the system, [61]. 

''Using the pseudopotential H;i3.4(l with g = ^^-a (/z = to/2) in the boson ground state Hartree-Fock 
mean-field equations H3.2.11|l leads directly to the Gross-Pitaevskii equation: 



Huang [73]: 




12 



(3.3.4) 
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(3.3.5) 



where n is the chemical potential and 0o is the (regular) single-particle wave function. 
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action which makes Vmf huge regardless of the other details of the potential. Remarkably, 
this leads to the counterintuitive conclusion that using a realistic two-body potential in the 
Hartree-Fock equations yields a much poorer result than using a 5-function potential with 
the same asymptotic scattering properties ^. This means, that it is not only convenient 
to make the shape-independent approximation in the Hartree-Fock approach but actually 
essential in the case of hard-core potentials in order to obtain quantitatively correct results. 
At the same time it is important to stress that the pseudopotential only works within the 
independent particle approximation, that is with the Hartree wave function, and should not 
be applied in an exact solution (see e.g. [63]). 

3.3.3 The Born approximation 

The Born series is the perturbation expansion of the scattering wave function or equivalently 
the scattering amplitude in powers of the interaction potential. It is interesting in this mean- 
field context because the condition of the pseudopotential to neglect the distortion of the 
(incoming) wave function in the region of the two-body potential, is precisely the same 
requirement that makes the Born approximation scheme valid in scattering theory, [3]. In 
particular the first term of the Born series follows directly from the assumption that the 
initial wave function is an undistorted plane wave, that is 

^L^^N ~ ^L^^lm = (27r)-3/V'^-^ (3.3.6) 

Plugging this into ()2.2.5|) gives the first-order Born scattering amplitude 

f{k', k)\B, = --^ / dxe~^^'^'~''^-^V{x) (3.3.7) 

which in the low energy limit, — > 0, yields the Born approximation scattering length 

as -/(O, 0) 1^, = ^ I dxV{x) = ^l^ drr'V{r) (3.3.8) 

where the last equality holds for central potentials. In the specific case of the zero-range 
pseudopotential ()3.3.4j) it follows that as = a, that is, the first-order Born approximation 
to the scattering length actually coincides with the s-wave scattering length. This fortu- 
nate property greatly simplifies the treatment of the N-body problem when one models a 
sufficiently weak interaction with the pseudopotential ^ [62,76]. 

It is also worth noting, that the value of is in general very different from the s- 
wave scattering length, a, for all potentials other than the pseudopotential. The numerical 
relation between a and as is illustrated in figure ITT] (^) for the cases of a Gaussian potential, 
a square well potential and the contact pseudopotential. When the former two potentials 

^In [62], B. D. Esry illustrates this by comparing Hartree-Fock calculations for the pseudopotential and 
for the realistic Morse potential with the exact hyperspherical result for three atoms in a harmonic trap. 

^In a more general treatment the Born scattering amplitude is replaced by the two particle T -matrix 
element which holds regardless of the interaction potential strength. In the low energy s-wave scattering 
case the T -matrix element is proportional to a, giving the same results as the Born approximation [6]. 

^Data for the Gaussian potential is numerically calculated (in atomic units) from subsequent runs: 
scatlen -VO Vp, where — 10~^ ^ Vq < 10~^ with steps of 10~^. For the square well potential one has 
analytically [2]: a — b{l — tan(7)/7), where 7 — y^2iib'^\Vo\/h'^. 
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Fig. 3.1: Scattering length, a, as a function of the Born approximation to the scattering length, a^, 
in units of 6 = 18.9ao, for a Gaussian two-body potential V{x) = Voe~^ {qb = ^/TTnb^Vo/2h'^) , 
a square well two-body potential V{x) = Vq, x < 6, V{x) = 0, x > h {ub = 2fj,b^Vo/3h^) and the 
contact pseudopotential V{x) = A'Kh^a5{x)/m {ub = a). 

are attractive (a^ < 0) the emergence of bound states and corresponding energy resonances 
are clearly visible. The similar overall behavior of the curves for these cases reflects the 
shape independence at low scattering energies. Obviously, the Born scattering length is a 
good approximation only when the interaction is very weak and a <^ b. 

3.3.4 Validity range of the pseudopotential approximation 

Since the condition for the application of the pseudopotential and the first Born approxi- 
mation are quite similar it is interesting to consider the validity range of the latter in more 
detail. One well known range where the Born Approximation describes the scattering am- 
plitude nicely is the high-energy regime, where the energy of the incoming particle is much 
greater than the energy scale of the scattering potential, [2]. This is not relevant for the 
low energy pseudopotential approximation described here. However, the requirement that 
the incoming wave function is not significantly altered in the region of the potential, or 
cquivalently that the second term in the Born perturbation expansion is very small, can be 
related via the Lippmann-Schwinger equation to the condition ( [3], p. 388) 

< 1 (3.3.9) 

For central potentials and at low collision energies {k 0, e'*^^ — > 1) this gives 

< 1 (3.3.10) 

which is obviously satisfied for a sufficiently weak scattering potential, V{x). The 5-function 
in the pseudopotential fulfills the condition by definition (0 <^ 1). In the case of a Gaussian 
two-body potential, V{r) = Vqc"'^ , the simple integral evaluates to 

lib'^\Vo\/h'^ < 1 (3.3.11) 

which is also the expression obtained for the square-well potential. This requirement may 
be compared with the condition for the potentials to develop a two-body bound state, that 
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is a solution to ()2.2.3p with E < 0, where E = —h'^/2fia'^ ( [4], p. 57). For the square-box 
potential the bound-state condition is /i6^|Vo|/^^ > vr^/S ~ 1.2 (by inserting the analytic 
expression for the scattering length, see caption of fig. 13. This is quite the opposite of 
condition (j3.3.1H) . A similar analytic expression for the Gaussian potential is not available, 
but in the case where b = 11.65 a.u., fi = 1.58 x 10^/2 a.u. (^^Rb), it is possible to determine 
numerically (by using the scatlen program) that the first two-body bound state occurs at 
Vo = —1.25 X 10"'^ a.u.. Plugging this into ()3.3.11|) gives fib'^\Vo\/h'^ ^ 1.34, which also 
clearly violates the condition. In other words, if the potential is strong enough to develop a 
bound state, the Born approximation and the pseudopotential approach will probably give 
misleading results. The conclusion is then, that a weak (attractive) interaction in the current 
mean-field theory, is one that is far from supporting a bound state. This observation is also 
visible in the numerical results presented in chapter El (see e.g. fig. 15. 8|) . 

3.4 Explicitly correlated description 

The key point of the preceding section is that the pseudopotential can be used under cer- 
tain conditions as a mean-field effective interaction and without the necessity of calculating 
detailed short range correlations. The conditions where found to be satisfied at low energies 
by weakly interacting dilute systems, where the particles are mostly far away from each 
other and correlations in head-to-head collisions are expected to be negligible. However, the 
importance of correlations must increase with the density of the system and the strength 
of the interaction, and at some point the mean-field approach becomes inadequate. Going 
beyond mean-field theory is only possible with the explicit inclusion of correlation effects in 
the wave function, that is, an appropriate choice of the correlation factor F in fl3.1.1|) . The 
following section describes a simple method for constructing F, in a way very similar to 
the discussion of translationally invariant clusters in coordinate space given by Bishop et al 
in [54]. 

3.4.1 Two-body correlations 

As the first step towards a systematic approach to the exact correlated ground state wave 
function one can consider a correlation factor containing only two-body correlations, e.g. 

N 

F2{r,,...,rM)=l[f2{nj) (3.4.1) 

where /2 is a properly chosen pair correlation function depending only on the interparticle 
distance ^ . This is the widely used Jastrow ansatz [24]. As stated at the beginning of this 
section the function /2 should go to unity, i.e. to the mean-field limit, at large separations 
manifesting the absence of correlations when the particles are far away from each other. 
At short distances the correlation function is expected to deviate from unity and writing 
f2{j'ij) = 1 + C2{rij), where C2{rij) represents the short range deviation, yields 

N N ^ N 

F2(ri,...,r^) = J](l + C2(r,,)) = 1 + ^ C2(ri,) + - ^ C2{r,,)c2{ru) + ■ ■ ■ (3.4.2) 

KT i<j ' {i<j)^(k<l) 

proper correlation function has to satisfy certain requirements (e.g. approach unity at large particle 
separation). See the complete list on page 62 in [8]. 
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where the indices of the interparticle coordinates, appearing in the summed products, at 
all orders, never overlap. The second term in this expansion corresponds to the effect of 
pair correlation while the third term induces separate correlations between two independent 
pairs of particles (clusters) and so forth. For a sufficiently dilute system it is unlikely that 
two or more independent pairs simultaneously are close in space and the expansion can be 
truncated after the first two terms, giving 

N N N 

F,{r^, . . . ,r^) ^ 1 + $^C2(r,,) = ^ ( _ + C2(r,,)) = ^^^^(r,,) (3.4.3) 

i<j i<j ^ i<j 

where C2{rij) is the simple redefinition of C2{rij) that absorb the factor one in the expansion. 



3.4.2 Three-body and higher-order correlations 

The most obvious improvement of the two-body correlation factor in ()3.4.3|1 is to include 
the next term in the expansion ()3.4.2|) . This is easily done by replacing C2(^jj) with a more 
general function, C^\rij,rki), depending on two interparticle distances, and realizing that 
it is possible to absorb all lower-order terms within this form, thus giving -^2(^1; • • • ^^n) ~ 
Yl!i<j^k<i^'2\^i3i'^ki)- In most cases, however, the corresponding improvement is small and 
the introduction of three-body correlations is much better (see [54], table 2). In particular, 
extending the Jastrow formulation to include three-body correlations, leads to 

N N N 

F3(ri, . . . ,rAr) = JJ/2(rij) JJ hijrij.rik.rjk) ^ ^ C3{rij,rik,rjk) (3.4.4) 

i<j i<j<k=l i<j<k=l 

where /s is a proper triplet correlation function and the freedom in choosing the functional 
form of C3{rij,rik,rjk) has been utilized to absorb all lower-order (cluster) terms. Following 
the same ideas, a strait forward generalization of the Jastrow approach (as done by Feenberg 
[77]) to include all higher-order correlations in the wave function gives 

FAr(ri, . . . , r^v) ^ 5 CAr(ri2, ris, . . . , r(7v_i)iv) (3.4.5) 

where the symmetrization operator, S, ensures that the correlation factor does not influence 
the exchange symmetry of the wave function, and the functional form C7v(?"i2, ''^is, • • • , ''^(Af-i)Af) 
is assumed to completely describe all the correlations of the N-body system. It should be 
noted, that in order to make a calculation manageable in practice it is necessary to further 
expand the unknown correlation functions, C„, in a set of simple functions (for example 
Gaussians, as described in detail in section 1^^ . 

3.4.3 Vahdity range of the Jastrow- Feenberg description 

Several points concerning the validity of the correlation description above are important: 

• Firstly, it is apparent that the application is limited to homogeneous and isotropic 
systems, that is, /i(Ti) = 1 and f2{ri,rj) = f2{rij), etc. The translational invariance 
resulting from this is essential to avoid problems with the center-of-mass motion [54]. 
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• Secondly, the correlations do not depend on internal quantum numbers such as spin. 
This is inappropriate in cases where the interactions are state dependent like in nuclear 
physics. Unfortunately, including state dependence in the Jastrow-type correlation 
functions, /j, turns them effectively into non-commuting operators demanding further 
symmetrization of the product form of F. This considerably complicates the formal- 
ism and is not considered here (for details see e.g. the FHNC single-operator-chain 
(FHNC/SOC) method [78] or CBF theory [8]). 

• Thirdly, the lack of momentum-dependency in the Jastrow-Feenberg ansatz makes 
it questionable when dealing with the ground state of Fermi systems, since there is 
no information about the specific location of a given particle within the Fermi sea. 
Such a treatment is perhaps acceptable for "integrated" quantities like the energy, 
but it is not at all clear whether it works for physical properties, like the specific 
heat, depending primarily on the "active" particles close to the Fermi surface [5]. To 
examine this question, it is again necessary to go beyond the Jastrow-correlated wave 
function (see [8], chap. 7). 

• Finally, returning to the discussion central in the pseudopotential approximation above, 
the particular choice of a wave function parametrization always corresponds to a re- 
striction of the full Hilbert space solution. While this restriction is quite severe in 
the mean-field Hartree wave function (with F = 1), leading to failure in combination 
with the exact interaction potential, it is much less pronounced with the Jastrow-type 
correlation factor. Still, the truncated factors, 6*2(^4^) and Cs{rij,rik,rjk), are clearly 
not able to take hard-core repulsion into account, since this requires that all pairs 
are simultaneously correlated. Only C]\f{ri2, ri3, . . . , r(^N-i)N) has this feature and thus 
seems to be valid with hard-core potentials. But whether a solution based on a realistic 
or, as in this work, a Gaussian interaction, with the given inclusion of two-, three- or 
N-body correlations, will reproduce reasonable results is, however, not obvious. The 
most convincing argument is, that many, if not most, of the key methods in modern 
quantum many-body theory are based on the Jastrow-Feenberg approach or similar 
ideas, and they reach such a high level of accuracy, also when limited to pair or triplet 
correlations, that it has been debated, although recently disproved (see [55]), if the 
ansatz could be generally exact. This obvious success story continues with the results 
presented in chapter El 
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Chapter 4 

The Stochastic Variational Method 



The SVM has been developed through the search for precise solutions of nuclear few-body 
problems [1]. In this chapter it is shown how to employ the method to atoms and N-body 
systems of trapped bosons. The main aim is to develop the SVM formalism to a point where 
one can systematically include the effects of two- and higher-order correlations in a way 
which is both intuitive and computationally efficient. Subsequent sections treat the trial 
and error selection procedure, the explicitly correlated basis functions and the details of how 
to symmetrize the trial function in practice. The three-body system constituting the Helium 
atom is used to benchmark the method towards treating the more intricate case of BECs. 

4.1 Stochastic trial and error procedure 

The variational foundation of the time-independent Schrodinger equation presented in chap- 
ter |21 provides a solid and arbitrarily improvable framework for the solution of diverse bound- 
state problems. A key point is that the quality of a variational calculation crucially depends 
on the trial function, \1/ = YliLi (^ii^i: consequently on the choice of the basis functions, 
{ipi, . . . j'iPk}- Assuming that each basis function depend on a set of nonlinear parameters, 
ipi = ilj{ai \ . . . ,ap^),i = 1..K, the SVM attempts to set up the most appropriate basis 
by a stepwise strategy: One generates a would-be basis function by choosing the nonlinear 
parameters randomly, judges it's utility by the energy gained when including it in the basis, 
and either keeps or discards it. In turn, each parameter is then changed (still randomly) in 
the search for additional improvement. One repeats this "trial and error" procedure until the 
basis found leads to convergence and no further energy is gained. The control flow structure 
for this method is shown in figure 14.11 (next page) and corresponds to a detailed version of 
the optimization fiow diagram displayed on the left of figure 12.31 This selection strategy has 
several advantages, where the most important are: 

• The optimization iteration (the fiow within the dotted box) is clearly separated from 
the computationally demanding solving of the generalized eigenvalue problem. This 
means that the nonlinear parameters can be improved repeatedly without the need of 
diagonalizing a i^'-dimensional matrix. 

• Due to the stepwise optimization procedure, a relative small number of matrix elements 
have to be calculated to test a new basis function candidate and the corresponding 
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K = 0, Eo = oo j 




Stochasticly select t trial values {cii, . . . , at} for the i'th parameter 






Calculate the t ground state energies {£'1, . . . , Et} that correspond 
to af~^^ = ai, . . . ,at in the trial wave function * = J2k^^ ^(a^+^) 






Locate the index min of the lowest energy from {Ei, ... ,Et} 


No^^ 


<-EoJ> 
Yes 



Choose af^^ = OLmin and set £^0 = E^, 



Increase K and solve the generalized eigenvalue problem He = eSc 




Optimized wave function is ^ = X^j^"*""" (l)k giving energy Ep ^ 
Fig. 4.1: The control flow diagram for the Stochastic Variational Method. 
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ground state energy is easily determined by finding the lowest root of a simple equation 
(see below). 

• The variational principle (j2.5.1|) ensures that the energy of the i^T-dimensional basis 
is always lower than that of an {K — l)-dimensional one. The procedure is therefore 
guaranteed to lead to a better and better upper bound of the ground state energy. 

Even though it is rarely the case, one still has to make sure that the solution is not on the 
plateau of some local minima. This is most easily done by confirming that independent 
calculations starting from different first basis states (i.e. different random seeds) lead to 
practically the same solution. 



4.1.1 Gram- Schmidt diagonalization 



For the stochastic optimization to be practical, it is essential that the ground state energy 
corresponding to a trial function candidate is evaluated with minimal computational effort. 
Otherwise it is simply not possible to test enough candidates to cover a reasonable parameter 
spread. A full diagonalization performed by solving the general eigenvalue problem is out 
of the question. Fortunately this can be avoided if only one basis function, let's say V'x+i, 
is changed or added at a time and the eigenvalue problem. He = eSc, for the other basis 



functions, ijji = ijj{ct 



1..K, has been solved. The idea is to evaluate the eigenvalues 



of the {K + l)-dimensional problem in a basis of orthonormal functions. Obviously, the K- 
dimensional solution has produced eigenvalues €1,62, ■ ■ ■ ,eK and corresponding eigenvectors 
. , c^^-* satisfying Sc = 1, and can be written in standard diagonal form 



c(i),c(2) 



Ai 
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eKj 



do 



ydK/ 



d2 
^dxj 



(4.1.1) 



k}, where 0i = Yl^icTi^j^i = ^--K. The 



in a basis of orthonormal functions, 02, • • • , s^a j , wi^v^xv. — >-j rj, 

(i^ + l)-dimensional solution can then be obtained by first applying Gram-Schmidt's method 
to construct (pK+i from ipK+i so that it is orthogonal to all 0i, 02, ••, 0_ft:, i-e. [9] 
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and then solving for the eigenvalues of 
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(4.1.2) 



(4.1.3) 



where hj = {(j)j\H\(f)K+i) and h* is the complex conjugate of hj. For this to work, the 
candidate ipK+i has to be linearly independent of the previous basis functions as required 
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by the Gram-Schmidt orthogonahzation ^. The characteristic equation of ()4.1.3|) is 



det{H-XI) 



ei - A 
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h2 



■ ■ ■ ex — X hx 

hi h2 ■ ■ ■ h*j^ hx+i — A 

which, assuming that all hi are nonzero, has a straightforward reduction 



(4.1.4) 



det{H-XI) = (-1)^+1/1* 
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(4.1.5) 



Thus, when (ej — A) 7^ (-v^ /ij 7^ 0), the + 1 eigenvalues Ai, . . . , Aj^+i are obtaining by 
simply finding the roots of the function 



K 



\hi 



(4.1.6) 



as graphically exemplified in figure ^21 The variational theorem ensures that ei < Ai < £2 < 
X2 < ■ ■ ■ ex < Ax+i, assuming both ei, . . . , e/^ and Ai, . . . , \k+i are arranged in increasing 
order. This is also helpful in setting up intervals for a root- finding algoritm (see app. ID.5|) . 

4.1.2 Refining process 

At any particular time, only one parameter out of the possible large number of nonlinear 
parameters in every of the K basis functions can be considered optimal, since the optimiza- 
tion procedure is applied consecutively, element by element, rather than simultaneously. It 
is then reasonable to expect that at least some of the previously added basis functions are 
no longer optimal or even needed. Especially when the basis functions are nonorthogonal to 
each other this might be the case, since, even though none of them are really indispensable, 
any of them can be omitted or changed because some others will compensate for the loss. To 
help shake off these flaws a so-called refining process is introduced. After having successfully 
found Kmax basis functions, one can further improve the energy without increasing the basis 
dimension. This is done by iterating through the current functions still optimizing their 
parameters in the spirit of the trial and error algoritm. After a few of these refining runs 
all of the basis functions play an active role again, and depending on the value of K, this 
process often improve the result considerably. 

^In practice, this must be explicitly verified during the stochastic selection procedure. 
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Fig. 4.2: Graphic illustration of the characteristic polynomial D{X) = X^j^j^ j^-^^ — X+hK+i, where 
the roots Ai , . . . , Xk+i determine the eigenvalues of the eigenvalue problem (|4.1.3|) . The example 
here corresponds to the adding of a candidate function ips to an existing basis {ipi, ■ ■ ■ , tpj} in the 
calculation of °°He (see section . There are four additional roots (A > 0) not shown here. 



4.2 Explicitly correlated basis functions 

From the outset, the SVM works well with any basis function that leads to analytical closed 
form evaluation of the required integrals of the Hamiltonian and overlap matrices. In ad- 
dition, the stepwise optimization of the variational parameters allows the efficient handling 
of a relatively large set of nonlinear parameters, . . . , ctp'*}, per basis function. This is 
important when trying to incorporate flexibility in the basis functions to treat complicated 
correlation effects. In the following it is described how to use explicitly correlated basis func- 
tions of the kind introduced in section 13.41 with the SVM. This treatment is applicable to 
both few-body and many-body systems, all though in practice, the full correlated description 
is only feasible for <~ 5. 

As discussed in chapter El the correlation description adopted in this thesis is based on 
the Jastrow-type trial function form, \1/ = -F(ri2, ris, . . . , r(7v-i)Ar)$(T'i, ^2, • • • fN), where F 
is the correlation factor and $ is the mean-field model state. Moreover, the factor F was 
approximated by symmetrized correlation functions, i.e. F ^ 5Cn(ri2, ris, . . . , r(„_i)„) in 
the case where up to n-body correlations are considered. To apply this with the SVM, it is 
necessary to expand the C^s in a mathematically complete set of functions where each term 
is simple enough to give analytical expressions for the matrix elements. Both Varga [1] and 
Wilson [22] argue that the only set of functions which meets such requirements for A^-body 
systems is the so-called contracted Gaussian basis ^ (i.e. Gaussians with different widths). 
For example, in the case of pair-correlation, this leads to the simple expansion 

C2(ri2) = J2 Cfcexp( - -ftfcr^^) (4.2.1) 

k 

where it is indicated that in practice the sum must be truncated at some finite level. In 



^Corresponding to the I = case of the nodeless harmonic-oscillator basis. 
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general, the expansion of the n'th order correlation function becomes 

C„(ri2, ri3, . . . , r(„„i)„) = ^ Cfc5 exp - - ^ "!^^^^^^,•) (4.2.2) 

k i<j=l 

where the symmetrization operator, S = P, includes perturbation terms for the 

first n particles only. There are many, possibly an infinite number of expansion sets, which 
approximate a given function by Gaussians equally well ^ . This makes the Gaussians an 
appropriate basis for stochastic optimization. 



4.2.1 Correlated Gaussian basis 

If expressed using the Gaussian expansion ()4.2.2|) . the N-body variational trial function can 
be written in the desired form, \E' = Ckipk-, where the basis functions are given by 

1 ^ 

7/>fe = $(ri,r2,...r^)5exp(-- ^Vj) (4.2.3) 

i<j=i 

and the set {ipk} is both non-orthogonal and over-complete (i.e. satisfying the requirement 
that results obtained by a systematic increase of the number of basis functions, will converge 
to the exact eigenvalues). These functions correspond to symmetrized sums of the explicitly 
correlated Gaussians originally suggested in 1960 independently by Boys [25] and Singer [26]. 
Over the years, they have been demonstrated to be an excellent basis for high accuracy 
variational calculations of few-body problems [1,27,32-34]. 

To utilize the correlated Gaussian basis in a translationally independent N-body solution 
it is necessary to write eq. ()4.2.3|1 in terms of the independent Jacobi coordinates defined 
by the matrix ()2.3.5|1 . Inverting the linear transformation ()2.3.3j) . one has 

ri = {u^^fx, and = Vi - = {u^'^^f x, (4.2.4) 

where u^^^^ = u^^^ —u^^^ and the vectors u*^*^ have components uf' = {U^^)ik- The summation 
in the exponent of the correlated Gaussians can then be written 

N N N-1 N-1 

i<j=l «<i=l k=l 1=1 

= x'^Ax, (4.2.5) 



where 



N 



A= (4-2.6) 

i<j=i 

TtP = u^^^\u^^^Y = {{U-%k - {U~%km-%, - {U-%i). 

The matrices T*-*-'-' with i,j = 1, . . . , A^, and hence A, are symmetric [N — 1) x [N — 1) 
matrices. Assuming that the model state, $, can also be expressed in Jacobi coordinates as 

■^Heuristic discussions on the completeness and fast convergence of Gaussians can be found in C6.1 of [1], 
the appendix of [36] and in [38]. 
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$(ri, . . . , riv) — > ^int{xi, . . . , XN_i)^cm{xN), the correlated Gaussians in the form ()4.2.3p 
are equivalent to 

ipk = ^int{xi, . . . ,XN-i)^cm{xN)S exp - ^x'^A^'^^x^ (4.2.7) 

where the independent entries of A are related to the aij parameters via ()4.2.6|) and 

aij = -iU''AU)ij (z<j) (4.2.8) 

One should note, that a necessary condition for the correlated Gaussians in ()4.2.3p and 
()4.2.7|1 to be square integrable, and thus have a finite norm, is that the parameters are 
positive and A is positive definite (i.e. x^Ax > 0, for all nonzero vectors x G M^^-^)) [16]. 
This must be explicitly checked when trying to optimize ipk by "guessing" with the SVM. 

The success of the correlated Gaussian basis is mainly linked to the amazingly simple 
form and consequently fast computation of the resulting matrix elements (the expressions 
are evaluated in detail in appendix [HI). Moreover, the form ()4.2.3|1 . with explicit correlation 
parameters a^j, presents a natural physical interpretation. For a one-dimensional Gaussian 
function, e~^"^ , the position expectation value is (r) = 2/ y/mi. Hence, when using the form 
()4.2.3|) and aij as variational parameters, 1 / ^y<yij can be viewed as an average "distance" 
between particles [1]. This makes ()4.2.3|) advantageous when setting up initial values for a 
variational procedure or in random selection, since valid intervals for the particle distances 
can be estimated from physical intuition (bound or trapped particles are not expected to 
move far away from each other!). 

4.2.2 Correlated exponential basis (A^ = 3 only) 

The Gaussian expansion is not always economical in describing the asymptotic behavior of 
the wave function at large distances. Only for Gaussian-type interaction potentials, har- 
monic oscillator potentials and in a few other cases does the exact wave function have a 
Gaussian asymptotic dependence. In the case of Coulomb systems and a wide number of 
similar potentials, e.g. exponential and Yukawa-type, the wave function has an exponential 
asymptotic. Moreover, the Gaussian expansion does not give a correct value for some specific 
short-range quantities such as the Kato cusp condition (see [31]). It is therefore interesting 
to consider the SVM with a trial function based on a correlated exponential basis. 

An exponential basis is, however, not amenable to analytical evaluation of the matrix 
elements for a system of more than three particles [46]. Consequently, only the = 3 case 
will be considered here. The application of the exponential expansion, gives 

^max 

Csiru, ri3, ^23) = ^ CkS exp(-afc?^i2 - f^kTn - '^kr2z) (4.2.9) 

k 

for the triplet-correlation function. Written in terms of the interparticle distances, denoted 
by a;^ = (^^127 '^13, '^23) for notational convenience, the corresponding correlated exponential 
basis functions becomes 



i^k = ^int{xi,X2,X3)S exp{-a^^^'^x) 



(4.2.10) 
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where a^*^^^ = {ak, Pk,'yk)- The details of calculating the matrix elements for this basis can 
be found in appendix IC.31 During the S VM optimization procedure, the inverse parameters 
and 7^^, are advantageously selected from those intervals in which the average 
distances between particles is expected to vary 

4.3 Symmetrization 

As outlined in section l^.fi. li the variational wave function should be either totally symmetric 
(bosons) or totally antisymmetric (fermions) under the interchange of identical particles. In 
the description leading to the correlated basis functions ()4.2.7|) and ()4.2.10p . the model state, 
$j„t, is assigned the role of imposing the proper overall symmetry. The most convenient way 
to achieve this with the SVM is by operating on the basis functions with the proper sum of 
permutation operators as defined by V in ()2.6.2|) . The specific permutations included in the 
sum of V depends on which particles are identical. 

In practice, it is helpful to represent a permutation P = {pi,P2, ■ ■ ■ ,Pn) of particle 
indices by a matrix, Tp, having elements [1] 

(Tp),, = 5,p„ z,j = l,2,...,Ar (4.3.1) 

Then, in the case of N identical particles, there are A^! different Tp matrices of size N x N 
and a specific permutation of the single-particle coordinates, P : rt ^ Vp., is simply written 

Pr = Tpr (4.3.2) 

Applying transformation ()2.3.3|) . the corresponding permutation will induce a linear trans- 
formation of the relative (e.g. Jacobi) coordinates, given by 

Px = fpx, where fp = UTpU'^ (4.3.3) 

Since the center-of-mass coordinate is unchanged under coordinate permutation and hence 
can be ignored, the size of Tp is (A^ — 1) x (A^ — !)• This way, the effect of P operating on 
the explicitly correlated functions, ipk, in ()4.2.7|) and ()4.2.10j) can be reduced to the simple 
replacements 

P : A^''^ f^A^'^^fp and P : a^^^^ a^^^^fp (4.3.4) 

With this approach, one can operate on ip^ with the symmetrizer S = P, where 

all A^! permutations are included in the sum, in the cases where the basis function is re- 
quired to be symmetric. If the spatial function is to be antisymmetric one should use the 
antisymmetrizer A = '^p{~^yPy where p=0,l is the parity of permutation P. In both 
cases, however, this will produce A^!^ terms in a single matrix element calculation. This gen- 
eral number is readily reduced to A^! since the correlated basis functions are of the product 
form = Y[i<j (with (pij being either a Gaussian of an exponential) so that any 

symmetric operator O satisfies {Pijipk'lOlipk) = {ipk'lOlPijijjk)- All matrix elements for Vipk 
functions can then be written 

{VMO\Vipk) = -^(MOlY^apP^k) (4.3.5) 
■ p 

''Follows from the discussion at the end of section 14.2.11 since in the case of a normalized exponential 
function, n{a)e~^"'^ , the expectation value of the distance, r, is (r) ~ a~^, [7]. 
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where ap = 1 for bosons and ap = (—1)^ for fermions and the single sum is over all N\ 
permutations of identical particles for both V = S,A. Unfortunately, the matrix element 
()4.3.5|) is still an 0{N\) computation which is only tractable for few-body systems (A^ <~ 5). 
Further simplifications of ipk has to be assumed to handle many-body problems (see e.g. 
section I4.5.3|) . 



4.4 Few-body system: The Helium atom 

The nonrelativistic ground state energy of the Helium atom has been a benchmark test 
for three-body calculations since the pioneering work of E. A. Hylleraas [21], 75 years ago. 
Recently, this subject has attracted much attention [45,47-51] and significant progress has 
been made, with the accuracy of the energy now at 36 decimals [50]. As a brief illustrative 
example, the SVM is now applied to the (l^S) ground state of °°He. The numerical results 
is presented in the next chapter and used to test the implemented computer program and 
the rate of convergence. 

The Helium three-body system consists of two indistinguishable electrons (labels 1 and 2) 
and an a-nucleus (3). Neglecting relativistic effects, the two-body interaction is exclusively 

Coulombic, with y^^^-Va = — —. The Hamiltonian (12.3.101) written in relative 

coordinates x'^ = (r'i2, ''"13, ■'"23) is then 

Ha nia \ ^ (4-4.1) 
_1 ^ 

ma Ma / 

where is the mass of the a-nucleus and fia = (^i!^^ ) is the reduced mass. In the present 
calculation we use = 00 making /Xq, = 1, ^. Following the approach described in the 
previous sections the trial function for the Helium ground state can be constructed from 
exponential basis functions: 

K 

^ = ^ CfcVfe, V'fc = ^{Xoo exp{-akXi - (3kX2 - ik^z)} (4.4.2) 

k=l 

where A is the antisymmetrizer, xoo = "^{<^(l)/5(2) — /3(l)ct(2)} is the two-electron singlet 
spin function arising from the coupling of two spin-| particles [7] and ak,Pk and 7^ are 
nonlinear parameters. An angular part in the trial function is not necessary for the L = 
ground state calculation. Since Xoo is constant although antisymmetric under the interchange 
of the identical electrons, it can be omitted if the antisymmetrizer is changed to A —>■ 1 + Pu 
(= 1S12), where P12 denotes a simple exchange of labels 1 and 2. 

All the necessary matrix elements are evaluated in appendix lC3l Since the basis functions 
are chosen to be real, both the overlap and Hamiltonian matrices are symmetric and the 
secular equation. He = eSc, can be solved effectively by well-known linear algebra methods 
[15]. The lowest of the eigenvalues found, ei, will then be the approximate ground state 
energy of Helium. The quality of the result will depend on the specific values chosen for the 
nonlinear parameters, ak,Pk and 7^, and the size of the basis, K. 

^Alternatively one might use the exact value TOq = 7294.2618241. 



1-T - 1 2 2 
H = --V^AV^ + , 

2 xi X2 X3 



with A 
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4.5 Many-body system: Bose-Einstein Condensate 

The main goal of this thesis is to discuss correlations in three- and four-boson systems. In 
the following is described how to employ the SVM to a system of identical bosons trapped 
by an isotropic harmonic oscillator and interacting via two-body potentials Vij. Moreover, 
four different levels of correlation are explicitly allowed in the variational trial function, 
ranging from mean- field to the full A^-body correlated treatment (as derived in chapter E)). 

Expressed in terms of the Jacobi coordinates defined in ()2.3.4|) the many-body Hamilto- 
nian describing the internal motion of a trapped N-boson system can be written 

_ _ ^2 . N 

Hi, 



where m is the boson mass and cu is the trapping frequency and 

Hem = -^^l. + \Nmu^W^ (4.5.2) 

is the center-of-mass Hamiltonian. It is apparent, that Hem represents the standard form of 
the three-dimensional harmonic oscillator having ground state energy Ecmfl = [7]. The 
total BEC energy is Eq = Eint,o + Ecmfl, where Eintfl is to be calculated by the SVM. 

Numerical computation in harmonic oscillator units 

In numerical calculations with limited precision arithmetic the optimal average order of 
magnitude of the numbers handled is ~ 1. To meet this demand in the case of BEC 
problems it is convenient to abandon the atomic units and do the numerical computation in 
the harmonic oscillator units, given by 



eho = fuj, aho = \l — (4.5.3) 
V mui 

where eho is the unit energy and a^o the unit length. With all lengths (rj, b, etc.) in units of 
ttho and all energies {H, Hem, Vij, etc.) in units of eho, the BEC Hamiltonian (j4.5.H) becomes 

1 '^^^ 2 ^ 
H^nt = -2 E - + E (4.5.4) 

i=l i<j 

and the ground state energy of a noninteracting trapped gas is just Eq = ^Neho- This 
means, that for a reasonable number of particles < 10^, the magnitude of the results in the 
harmonic oscillator units are also reasonable (~ A^). However, care must be taken during 
evaluation of the matrix elements since they can reach much greater values and are the main 
source for loss of accuracy. For very large N and in some other cases additional rescaling is 
required ^. 



'A technique for scaling the magnitude of the overlap (ipk'lipk) is demonstrated in appendix ID. 21 
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4.5.1 Selecting the BEC ground state 

In the case where the bosons are interacting attractively, the eigenstate with lowest energy 
is not necessarily the BEC state we are looking for. If the scattering length is large enough 
and N > 2, the bosons may form "molecular-type" many-body bound states even when 
the boson-boson interaction potential is too shallow to support two-body bound states. 
These states could as well be characterized as condensed (N-body) states but they exists 
only at high densities making them unstable to recombination processes. In addition, such 
bound states do not have the distinctive BEC features (e.g. density profile) obtained in 
experiments [52] and it is therefore important to select the correct "gas-like" condensate 
state as the target of the SVM. 

The characteristic difference between the self-bound and the trapped condensate in the 
attractive boson system is their spatial extension. To illustrate this, it is convenient to 
examine the effective potential experienced by a boson as a function of the hyperradius, p, 
i.e. the average distance between the bosons in the trap (see definition (j2.3.6|) ). A rough 
sketch of the behavior (details can be found in [64,66]) for the = 10 case is outlined in 
figure 14.31 a) (solid line) and shows a global minimum at low p and a second minimum 
(with U{p) > 0) at hyperradii around the center of the trap {ptrap = &t a/3A^/2). The 
effective potential for the corresponding repulsive interaction (dashed line) has only the 
second minimum and is almost indistinguishable from the solid line in the bottom inset. 
This second minimum supports (quasi-stationary) states with the characteristic features of 
Bose-Einstein condensates and the lowest of the corresponding eigenstates is the BEC ground 
state of interest here. 

Considering the attractive boson system in detail, it becomes apparent that the barrier 
enclosing the BEC minimum gradually declines as the number of particles or the scattering 
length increases. This is demonstrated in graph b) and c) of fig. 14.31 The barrier completely 
disappears roughly when \a\N/bt >~ 0.67, as derived previously by many authors [6,66,75]. 
This is the well-known limit where the Gross-Pitaevskii mean-field theory breaks down. 
However, in the current method, the BEC eigenstate does not collapse when a is increased 
to where the barrier vanishes, but instead, transcends smoothly down the potential hill, to 
become a weakly bound so-called Efimov state ^. At the same time, increasing the scattering 
length also makes the first potential minimum deeper and, independently of the forming 
Efimov state, allows more and more (lower lying) molecular- type bound states to appear. 

Taking these conditions into account, the SVM has to be targeted to calculate the energy 
of the BEC state in a special way. From the outset, the trial and error procedure can be 
designed to minimize the variational energy, ej, of any given eigenstate, in accordance 
with the general solution ()2.4.1|) and the variational theorem fl2.5.6p . However, to be able 



^The model potential used for this graph is derived in [64] by use of the adiabatic hyperspherical expansion 
method and composed of terms for the external trap (~ p^), the generalized centrifugal barrier (~ p~^) and 
an interaction part from the angular equation (~ Mp)) s-s 

^(ri^|lf^+ '^"-yf ,45,5) 



where bt = 0,^0 = ^Jh/rauj is the trap length. 

^The many-body Efimov states are unavoidable for large scattering lengths and are located in the barrier- 
absent plateau region of fig. 14.31 c) (for a = ±00), far outside the range of the two-body interaction but 
before the confining wall of the trap [71]. 
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Fig. 4.3: The effective boson-boson potential (|4.5.5j) as a function of the hyperradii for a Gaussian 
two-body interaction, V{x) = Voe~^ , with b = 18.9 a.u. a) The case of a weak attractive and 
a weak repulsive interaction for = 10. The insets show the finer details of the barrier region 
and the trap center, b) The details of the barrier region when the scattering length is fixed at 
a = —17.6 a.u. and the number of particles is varied and c) when the scattering length is varied. 



to specify the number, i, of the BEC eigenstate, requires that one knows the exact number 
of lower lying self-bound states (including Efimov states) before the calculation. All though 
crude estimates for the number of bound state exists (see [64]) they are not nearly precise 
enough to be applicable. The best alternative, which has been chosen from the many different 
schemes tested during this work, is to have the SVM always minimize the lowest positive 
eigenvalue. This means, that the algoritm will determine, at runtime, what eigenstate of 
the current basis corresponds to the lowest positive eigenvalue, and add the particular basis 
candidate that minimizes this eigenvalue. The target state will then automaticly increase 
by one each time another bound state appears in the solution. This procedure is illustrated 
in practice in section IB. 2. 21 

4.5.2 Mean-field 

The single-particle wave functions, (f)Q{ri), for the ground state of non- interacting bosons 
trapped by a spherical symmetric external field, have a Gaussian form [6], i.e. 0o(^i) ~ 
exp{—rf/2bt), where bt = auo = \/h/muj is the trap length. Using relation ()2.3.6p . this leads 
to a Hartree mean-field wave function that can be expressed in terms of the hyperradius p: 

N N 

<^HF = Y[Mr^) ~ exp( -J2r^/2bt) = exp( - iVi?V2&t) exp ( - pV2&t) (4.5.6) 

1=1 i=l 
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where R is the center-of-mass coordinate. Since the hyperradius is defined by = Si<j '"ij' 
any exphcit dependence of p corresponds to a mean-field infiuence, that is, the same cor- 
relation between any pair of particles. As Boson-boson interactions modify the Gaussian 
shape of the non-interacting system, an appropriate mean-field trial function would be an 
expansion in Gaussians depending solely on p, e.g. 

*MF = Xl^fc^^' = exp(^ - -iVa^'^V^) (4.5.7) 

k=l 

where the center-of-mass dependence is explicitly removed. Transforming to Jacobi coor- 
dinates with = ^^7^ x^, these basis functions become simple editions of the explicitly 
correlated Gaussians (i.e. ()4.2.7|) with A'-'^-* = a^^^I and l>j„i = = !)• It is straight for- 
ward to insert = (a^'^'^ -^-a^^'^)"^! together with identities ()C.1.9|) in the matrix element 
expressions for the correlated Gaussians given in appendix Kill and find ^ 

2^ \ 3(iV-l)/2 



/, iTr /3(Ar-l)mcu2-a('=')aW N{N-1) ^ , 

Hk'k = {^H^ntm = [ 2 Na(^') + NaO') + 2 ^ (V2)jg.-. (4.5.9) 

for the overlap and Hamiltonian mean-field matrix elements. 



4.5.3 Two-body correlations 

In a dilute system of interacting particles one can often assume that, at any given time, only 
two particles are close enough to each other to interact [65]. The rest of the particles are 
only "feeling" the mean- field. In such a situation, most often defined by n\a\'^ ^ 1, where n 
is the density [52], one should expect two-body correlations to be the dominant interparticle 
relationship and an appropriate trial function to describe this would be 

K 

^2B = J^Cfc^fc, i^k = 5 exp( - -(/3(^-) - a«)r22)exp( - -Na^'^p^^ (4.5.10) 

k=l 

where all pairs have the same mean-field correlation parameter, a^''\ except one pair (particle 
1 and 2 in the first term of S) that are correlated by The symmetrization makes sure 
that all separate pairs are taken into account in the same fashion. This ipk is already in 
the favorable form of a sum of the explicitly correlated Gaussian basis functions ()4.2.3|) . In 
addition, due to the fact that they are so simple, it is possible to derive analytical expressions 
for the matrix elements of Hint, that are independent of the number of particles in the system 
(i.e. with computational complexities 0{1)). The formulas and further details can be found 
in appendix IC2I 



4.5.4 Higher-order correlations 

While the dominant effect of interactions in dilute gases when n\a\^ -C 1 is due to two-body 
encounters, three- and higher-order correlations should become more and more important 

^Expressions for v{l/2) can be derived from HC.1.10|l - (jC.1.12|l . 
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with increasing density, n, or scattering length, a. The specific SVM trial function that 
includes the up to m-body correlations (m < N) of these cases, can be written 

*mB = J]c,^fe, i^k = 5 exp( - - J]a(JVj)exp( - -Na^'^p') (4.5.11) 

k=l i<j 

where ipk is a symmetrized sum of explicitly correlated Gaussians. In this expression, the 
individual correlations of the (m+ l)m/2 particle pairs is represented by an equal number of 
nonlinear parameters a\j\ Unfortunately, a full correlated treatment with m — N requires 

A^! different term in S and limits it's usability to EEC's of only a few atoms (maximum of 
N — 5 in this work). 



Chapter 5 
Numerical results 



This chapter presents the numerical results obtained by applying the SVM to the N-body 
systems introduced in sections 14.41 and 14.51 Starting with the well explored °°He atom, the 
implemented computer program is first thoroughly tested and bench marked. It is illus- 
trated that the convergence rate is fastest when the basis functions represent the asymptotic 
behavior of the exact wave function well. The main calculations will subsequently treat 
BEC systems of sizes = 3,4 and 10, where the bosons are interacting attractively over a 
wide range of scattering lengths (e.g. — oo < a < 0). Detailed graphs of the lowest energy 
levels are presented and it is shown how to distinguish between "gas-type" and "molecular- 
type" eigenstates. In addition, each individual calculation is repeated three times with the 
SVM trial function including different degrees of correlation [mean-field, two-body and full 
correlations) giving a clear indication of the importance and effect of correlations in such 
systems. 

5.1 Method test: ^He 

The strengths and weaknesses of the SVM can be rigorously investigated by considering the 
°°He three-body system. For a definite basis size there is a total of 3K nonlinear parameters 
(see fl4.4.2|) ) which, ideally, have to be optimized. Postponing optimization for a moment by 
considering completely random parameters reveals the result of expanding the function space 
Vk by functions that are far from optimal. The graph on the left of figure 15.11 shows the 
energy convergence of Helium as the basis size is increased from 1 to 50 by adding exponential 
basic functions given in ()4.4.2|1 . where the inverse parameters, a^^,P^^ and 7^^, are selected 
randomly in the intervals [0,4], [0,2] and [0,2] respectively ^. The three curves correspond 
to three different random seeds. The graph on the right shows the same convergence in the 
case of a Gaussian- type-basis ^. 

It is apparent from the convergence shown in fig. 15.11 that the crude adding of linearly 
independent basis functions is actually an effective way to reach the accurate ground state 
energy. This is the case for both types of basis functions. The corresponding stepwise 
construction of the appropriate wave function is shown in figure where the coordinate 
s is the length of the vector from the center of mass of the two electrons to the a nucleous 

"'^ Because the average distance between the electrons is expected to be twice the average distance between 
an electron and the nucleus, the interval for a'^^, corresponding to ri2, is set to twice that of /J^T^ and 7^7^- 
^i.e. <f>fc = (1 + Pi2)exp(— ittfexf — \(}kx\ — \'ykX%)- See section 11.2 .11 for further details. 
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Fig. 5.1: The energy convergence of the °°He system as a function of the basis size with (left) 
random exponential basis functions and (right) random Gaussian basis functions. 



and it is assumed that -L s for the sake of illustration. Clearly, if the SVM trial function 
is flexible enough, the variational principle ensures slow but sure convergence to the optimal 
representation as the basis size increases. 

Further optimization of the nonlinear parameters will improve the rate of convergence 
and limit a possibly excessive use of basis functions. In figure IHT^ the results of the SVM trial 
and error optimization is illustrated. The accuracy in correct decimal digits is displayed as a 
function of the basis size in the cases where 1,2, 10 and 100 random values are tested for each 
nonlinear parameter before adding the best trial. Again, the graph on the left corresponds 
to the exponential basis (j4.4.2|) where the inverse parameters, ^ and 7^^, are selected 

randomly, and the right graph corresponds to a Gaussian basis where the squared inverse 
parameters, and 7^^, are selected randomly. To avoid frequent linear dependency in 

the basis the random number intervals are doubled, i.e. [0, 8],[0, 4] and [0, 4] respectively, for 
all calculations. Since the exact wave function of Coulombic systems will have an exponential 
asymptotic behavior at large distances [46], the exponential basis produces much better 
results than the Gaussian basis. More importantly, however, this means that in the case of 
BEG calculations, as treated in the following section, the asymptotic is expected to have a 
Gaussian form, and therefore the Gaussian basis would be the best choice. 

From the curves it is obvious that optimization of the nonlinear parameters improve the 
accuracy of the results to a certain extend. The positive effect saturates, however, when the 
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a) K=1 , E=-1 .1 9309 a.u. b) K=5, E=-1 .94853 a.u. 




Fig. 5.2: Graphic illustration of the normalized wave function, ri2S^(ri2, s), of °°He, where = 
+ (ri2/2)^ = rig + (ri2/2)^, as calculated by the SVM with a Gaussian basis of increasing size 
reaching better and better energies. The three cases (a),(b) and (c) with basis sizes K = 1,5 and 
10, corresponds to a completely random selection, while case (d) is the result of an optimized trial 
and error selection producing the exact ground state energy, E = —2.90372. 

flexibility of the basis functions, which is very limited in this case, is completely exploited 
by testing many different random values for the parameters. Still, the variational theorem 
guarantees better results with every increase in the basis size. This illustrates the generic 
trade-off between high optimization and large basis size. On the one hand, focusing on 
flexible basis functions with many nonlinear parameters and high optimization costs, gives 
good results even for very low basis sizes ^. On the other hand, keeping the optimization cost 
to a minimum by using simple basis functions, allows a huge basis size and correspondingly 
precise results. One of the best values (24 decimal accuracy) of the nonrelativistic °°He 
energy has been achieved by V. I. Korobov [49] using K = 5200 simple exponential basis 
functions like ()4.4.2|1 with complex parameters selected pseudo-randomly from optimized 
intervals. However, as explained in the next section, this is not an effective approach for 



•^Thakkar and Koga [47] reach an impressive 15 decimal accuracy in the °°He ground state energy with 
a basis size of only K = 100. 
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Fig. 5.3: The accuracy in correct decimal digits (i.e. logio(-g^'^^^r^)) of the nonrelativistic ground 
state energy of °^He as a function of the basis size when the basis functions are optimized by 
testing 1, 2, 10 and 100 random values for each nonlinear parameter. The left figure corresponds 
to exponential basis functions and the right figure to Gaussian basis functions. 



more complicated systems, like the case of trapped bosons, because of the much wider 
random intervals needed. 

One may say, that achieving the high accuracy obtained in this section is redundant 
and has no physical meaning. However, the extraordinary precision is a consequence of the 
variational stability of the energy eigenvalue and does not necessarily reflect that the correct 
analytical structure of the wave function has been found ^ . The method above gives much 
poorer accuracy for the calculation of observables other than the energy, e.g. relativistic or 
QED corrections [30]. Obviously though, the results show the power of the SVM and of 
modern computers ^ and their ability to solve at least quantum three-body problems to any 
number of digits. 



^The exact wave function must satisfy the Kato cusp conditions [35,45] and includes e.g. logarithmic 

terms, which have ncgUgiblc effect on the value of the variational energy [30]. 

^AU the results presented in this section were computed on an (old) Athlon-650 CPU using 32-bit floating 
point arithmetic in less than one hour. 
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5.2 Bose-Einstein Condensate 

This section presents the numerical results obtained from applying the SVM to specific 
systems of attractively interacting bosons in the case of a spherical harmonic trap. The 
default physical parameters used in all the calculations are: 

Mass of ^^Rb 
Frequency of trap 
Lenght of trap (= ato) 
Har. Osc. energy 

Two — body interaction 
Range of potential 
Depth of potential 
S — wave scattering length 
Two — body bound states 

Note, that in the case of the Gaussian two-body interaction, it is the potential depth which is 
changed in order to vary the scattering length, while the potential range is fixed. Moreover, 
the potential depth is limited to values where no two-body bound states are supported. 

Four different combinations of trial wave functions and two-body interaction potentials 
have been considered and is denoted by the following names: 

• Hartree: Denotes a SVM calculation with the Hartree single-particle product in (j3.2.H) 
as the trial wave function and the Gaussian potential as the two-body interaction ^. 

• Mean-field: This corresponds to Hartree case but with the zero-range pseudopotential 
as the two-body interaction. 

• Two-body: A SVM calculation with a trial wave function given by ()4.5.im) that explic- 
itly includes pair correlation and a Gaussian two-body interaction. 

• Full: This name designates a SVM calculation that explicitly allows up to N-body cor- 
relations in the variational trial function (see ()4.5.1H) ). and has the Gaussian potential 
as the two-body interaction. 

In all of the above cases, the random value interval from which the nonlinear variational 
parameters are selected, is given by 

(aSf G [0.0001; 10] h (5.2.1) 

However, since this interval spans an impressive 5 orders of magnitude in the attempt to 
account for both molecular-type bound states and gas-type BEG states, the random value 
generator has to be specially designed to output an equal number of values at each order 
(e.g. as many parameters selected in the range [0.0001; 0.001] as in the range [1; 10]) ^. 

^This does not correspond to a genuine self-consistent Hartree-Fock calculation since the range and depth 
of the interaction potential are not variational parameters in the current approach. 

^The simplest way to do this, is by choosing a random number, v, from the interval [—4; 1] and then 
assign the nonlinear parameters as (a*^'^-')"^ = logio v. 



m = 86.91 amu (1.443 x 10"^^ Kg) 
u = 77.87 Hz 
bt = 23024 a.u. 

eho = 1.183 X lO"^'^ a.u. (5.160 x lO"^^ j) 
V{r) = Vbe""'/^' or V{r) = ATxh^a5{r)/m 
b = 11.65 a.u. 

Vo = [-1.248261 X 10"^; 0] a.u. 
a = [-10^; 0] a.u. 
iV;, = 
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Many previous calculations on BECs suggest that reasonably dilute condensates are well 
described by the s-wave scattering length alone [74,75]. In addition, the validity of the widely 
used Gross-Pitaevskii mean-field theory, which has been recently examined depends on the 
factor 0?. Therefore it is convenient, also in the current context, to describe the properties 
of trapped N-boson systems as a function of the s-wave scattering length, a. 

5.2.1 System with = 3, and — oo < a < 

The translationally invariant three-body problem for identical particles has only two degrees 
of freedom in coordinate space. This means, that the application of the SVM with basis 
functions having two nonlinear parameters is sufficient to include all correlations of three 
trapped particles. Consequently, the restriction to = 3 leads to a computationally simple 
problem and presents an opportunity to shed some initial light on the characteristics of 
BECs. 

Energy levels 

The overall behavior of the internal energy levels as a function of the s-wave scattering 
length, a, for three ^^Rb atoms in a spherical trap is shown in figure 15.41 Data points for 
the energies, Eq,Ei,E2 and £'3, corresponding to the four lowest eigenstates, \E'o,\I'i,\E'2 and 
\l/3, are plotted. The upper graph shows the energies that fall in the vicinity of zero and the 
lower graph displays the energies which are large and negative on a logarithmic axis. This 
rather complicate energy level structure is interpreted as follows. In the case where a is very 
close to zero, the energy levels correspond to the non-interaction boson gas result, given by 
En = {3N/2 + 2n)fvjj. Increasing the attraction slightly to where a ~ —0.002 6t, creates the 
condition for a molecular-type three-body bound state, and the lowest energy, Eq, "plunges" 
downwards. The second lowest energy, Ei, then takes the place of the lowest gas-like energy 
which is not noticeably affected by the forming molecular state. This sequence of events 
is repeated when the scattering length is further decreased to roughly a ~ —0.04 6^. From 
here on, down to where a ~ —100 it is only \l/2 and higher eigenstates that give energies 
above or close to zero. Moreover, these energies change dramatically in the region around 
a ~ —1 where E2 becomes negative and E^ drops by almost "ifiw. 

As a first conclusion, it is apparent that in the specific range —100 6^ < a < (with N^^ = 
0), the = 3 system has at most two strongly bound states {Eq and Ei for large negative a), 
which are interpreted as true molecular-type states. Secondly, there are additional unbound 
eigenstates, which seem independent of the molecular states, and, in the limit of weak 
interaction, are equidistantly spaced by 2hu!, corresponding to an ideal trapped gas. These 
states are interpreted as gas-type BEC states. The significant change that the gaseous energy 
levels undergo around a ~ —1 64 is linked to the disappearance of the barrier in the effective 
boson-boson potential (see figure . As the barrier vanishes the energy level "located" in 
the trap minimum of the potential is free to descend down in the global minimum. At some 
point, when the scattering length is much larger than the size of the trap, the energy level 
stabilizes. The corresponding bound state is considered to be a so-called Efimov state since 
it is weakly bound and has a large spatial extension (see below). 

Before ending this subsection it is worth noting that the gross behavior of the A^ = 3 
energy levels presented here, agree both quantitatively and qualitatively with the results 
obtained by D. Blume and Chris H. Greene in [75] with the adiabatic hyperspherical expan- 
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-100 -10 -1 -0.1 -0.01 

a/bt 

Fig. 5.4: Energy levels for the = 3 system of trapped bosons as a function of the scattering 
length. The data points show the energies of the ground state, ^'o, and the three lowest excited 
states, ^1-^3, as indicated. The upper figure details the energy levels of the the gas-type BEC 
states and the lower figure illustrates the "plunging" molecular- type energy levels. 

sion method. Consequently, one can confidently assume that the SVM also produces correct 
results when continuing into the (unchartered) many-body regime in the following. 

Characteristics of the BEC state 

It is of interest, also in the current context, to examine the defining features of the BEC 
eigenstates in more detail. In the simple = 3 case, one can illustrate the spatial dependence 
of the calculated wave function graphically as a function of the two degrees of freedom (e.g. 
ri2 and s, see fig. 15. 2p . This is done in figures 1375115 . 71 for three different scattering lengths. 
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Fig. 5.5: Graphic illustration of ri2S^Q,ri2S'^i and ri2S^'2, as a function of ri2 and s (see caption 
to fig. 15. 2() . where ^'o, and ^'2 are normalized and correspond to the three lowest eigenstates 
determined by the SVM for = 3 and o = —50 a.u. (= —0.0022 bt). These eigenstates have 
nonzero regions at ri2 ~ s ~ VNbt only, and are consequently interpreted as the ground, first and 
second excited BEC states. Notice the logarithmic scale on the ri2 and s axes. 
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Fig. 5.6: As above, but for scattering length a = —57 a.u. (= —0.0025 bt), where the first molecular- 
type bound state is close to appearing. At this scattering length, "ifQ and ^'i still correspond to 
the ground and first excited BEC states since almost all the amplitude is concentrated in the 
~ s ~ ^/Nbt region (N.B. the logarithmic axes). The ^'2 eigenstate, however, now resembles an 
unbound {E > 0) molecular- type state with the amplitude distributed at very small interparticle 
distances, ri2 ~ s ~< 10^^ bt- 




Fig. 5.7: As above, but for scattering length a = —100 a.u. (= —0.0043 bt). In this case, the lowest 
eigenstate, '^q, is a true molecular-type bound state almost unaffected by the external trap. The 
first and second excited eigenstates, ^'i and *I'2, correspond to the ground and first excited BEC 
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In addition, the corresponding energy and the result from calculating the root-mean- 
square radius, given by 

N K K N-1 

< rLs >= (^1 - R)Vm^) = {^\pVn\^) = J2Y1 E^^^i^'i^^^') (^.2.2) 

i=0 k=0 k'=0 i=0 

is indicated above each image. The expected and apparent conclusion from these graphs 
is that there are distinctive differences in the spatial distribution of the particles in the 
bound molecular-like eigenstates and the gas-like eigenstates. This significant difference in 
the spatial extension, or equally in the density, is then an effective measure to determine the 
type of a given eigenstate calculated by the SVM. In the following, it is therefore assumed 
that an eigenstate having Vrms >~ 10~^ h is a BEC type state. Look in the captions for 
further details of the interpretation. 

Correlations 

In reference to the above assumption, the reader might wonder about the wave functions in 
figure 15. 6[ where the amplitude is separated in two distinct peaks corresponding to different 
densities. Whether this mixture of wave function amplitude characteristic to both gas-like 
and molecular-like states is an actual physical fact or just a consequence of the stochastic 
method employed, is not clear from these calculations. It is quite possible that the peak at 
low interparticle distances in the first two illustrations of fig. 15.61 is a remnant of previously 
optimal basis configurations in the stepwise trial and error procedure. If this is the case, 
further rigorous optimization should slowly remove the high density peaks. In the discussion 
of correlations below, however, this is not really relevant as long as the root-mean-square 
distance gives a clear indication of the type of a given eigenstate. 

With the preliminary treatment completed it is now possible to turn the attention to 
the interparticle correlations of the three-boson system and concentrate on the effects they 
produce in a particular eigenstate. The most convenient way to illustrate these effects is by 
isolating the interaction energy contribution to the total energy, since, in accordance with 
the definition adopted in section ITTl this equals the correlation energy, Ecorr, apart from a 
sign. In general, for the ground state of N interacting bosons in a spherical symmetric trap, 
one has a simple relationship, given by E^orr = ^Nfvjj/2 — Etotau where the term 3Nhu!/2 
represents the energy of the non-interacting gas. 

Figure 15.81 shows the correlation energy as a function of the scattering length for the 
lowest BEC state of the = 3 boson system. Four different sets of data points are plotted, 
corresponding to the correlation levels included in the Hartree, Mean-field, two-body and full 
SVM calculations. The solid line, representing the full correlated treatment, is here regarded 
as the closest candidate to the "exact" correlation curve. As noted above, there is only two 
independent spatial coordinates in the = 3 case, and consequently, the two-body and the 
full correlated treatments should be equivalent. Evidently, this is also the case, in that the 
dashed line of the two-body correlated calculation is indistinguishable from the solid curve. 

The most apparent feature in fig. 15.81 is the obvious failure of the mean-field treatment 
to account correctly for the correlations in the system when the scattering length becomes 
large. This can be directly related to the break down of pseudopotential approximation 
as was predicted in section Hi. 3. 41 In the recent article [67], DuBois et al establishes the 
condition, (A^ — l)|a/6tp <~ 10~^, for the validity of the GP mean-field calculation of 
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Fig. 5.8: Correlation energy of the lowest BEC state, -BbecO; in t^ie A'' = 3 boson system, as 
a function of the scattering length for SVM calculations including different levels of correlation. 
Note that the data points produced by the full treatment and the two-body treatment are almost 
identical up to the precision considered here, and consequently indistinguishable in the graph. 



bosons in a trap. With = 3 this amounts to a >~ —0.08 bt in the current case. Looking 
at the curve for the mean-field SVM calculation, this rough limit agrees very well with the 
observed validity range. 

One last lesson may be learned from comparing the Hartree and the mean- field curves. 
Both of these treatments are based on the Hartree product trial function which is explic- 
itly uncorrelated and the only difference is in the adopted two-body interaction potential. 
However, the results are very distinct and the Hartree calculations using the finite Gaussian 
potential as the two-body interaction is clearly an unlucky choice which leads to a terrible 
treatment of correlations. The mean-field calculation, on the other hand, reveals that the 
application of an effective interaction as opposed to a realistic one, can be quite powerful in 
the attempt to include correlation effects. This might explain why the mean-field approaches 
have been exceptionally successful in treating dilute systems of bosons. 



5.2.2 System with iV = 4, and -oo < a < 

The system of four trapped bosons, in particular, has the features to be a great source of 
knowledge in a discussion of many-body correlations in BECs. With = 4 particles it is 
expected, that the interparticle correlations governing the dynamical motion will include 
both three- and four-body effects. At the same time, such a system is simple enough to 
facilitate numerous and accurate calculations of the supported eigenstates and corresponding 
energies. In the present application of the SVM, these features makes the = 4 system 
unique, since a subsequent increase in the number of particles makes the full calculations 
intractable. 
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Fig. 5.9: Left: Energy levels of the two lowest BEC states for the = 4 system of trapped bosons 
as a function of the scattering length. Right: The corresponding number of lower lying (with 
respect to the BEC state) molecular- type states, both in the case of the full SVM calculation and 
the two-body SVM calculation. 



Energy levels 

Energies corresponding to the two lowest BEC states, "^becfi and "^bec,!, of the = 4 system, 
have been calculated with the SVM in the two-body description. The resulting levels are 
displayed in left part figure 15.91 In the same graph, the energy curve obtained with the full 
SVM calculation is also plotted, however, only for the lowest BEC state and for a limited 
range of the scattering length (—2 6* < a < 0). Apparently, the characteristic behavior of the 
energy levels is very similar to the = 3 case in figure As is the related interpretation. 
One slight change can be observed by taking a closer look at the range where Ef,ec,o falls 
off significantly. This now happens around a ~ —0.3 bt, and indicates that the barrier in 
the effective boson-boson potential vanishes at a smaller scattering length than in the three- 
boson case. This comes as no surprise since the disappearance of the barrier is known occur 
when \a\ ^ Q.67bt/N, [66]. 

On the right hand side of figure 15.91 a step-graph shows the number of molecular-type 
bound states determined by SVM calculation as a function of the scattering length and 
corresponding to the data points in the left figure. Contrary to the = 3 system, there 
are now a large number of lower lying bound states to take into account. However, this 
depends heavily on the correlations allowed in the trial wave function. While the two-body 
treatment finds at most five "plunging" energy levels, the full calculation determines up 
to 28 in the reduced scattering length range considered, and even more appear for larger 
attraction. As is pointed out later, this is the reason why the full calculation has been 
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limited to scattering lengths above —2 ht- Moreover, the significant difference in the number 
molecular-type states supported at particular levels of correlations was very much expected, 
since the density for these states is so high, that higher-order effects should play an important 
role in the interparticle dynamics. The data presented here confirms this prediction. 

In light of the main topic of this thesis, the most interesting feature visible in the = 4 
energy level figure, is the revealing separation between the two-body curve and the full 
curve in the region where they are both calculated. Clearly, the full treatment results in 
a significantly lower total energy than the two-body treatment when the scattering length 
becomes lower than ~ —0.2 bt- In other words, a description including all higher-order 
correlations, produces a distinctively better variational upper bound to the BEC ground 
state energy, than one including only two-body correlations. This discovery is very important 
and will be further investigated below. 

Density profile 

One of the defining properties of a BEC is the characteristic density profile of the condensed 
gas. To further support the understanding of the lowest gas-like state from the SVM cal- 
culations as a true BEC state, it is therefore convenient to consider the one-body density 
function, which is given by [54] 

N K K N-1 

n{r) = (v&l 5^5(ri - - r)|M/) = J] J] ^(^fc|5((t.«)^x, - r)|^,,) (5.2.3) 

i=0 k=0 k'=0 i=0 

where = {u^^^)'^x and R = x^^ was used ^ . Figure l3'.10l illustrates this function for the low- 
est gas-type eigenstate determined by the SVM in three specific cases of increasingly negative 
scattering length. Obviously, the solid line, calculated for very low a = —0.00434 bt = —100 
a.u., lies very close to the form of the analytically available Gaussian shape (n(r) oc e"'' ) of 
the non-interacting ideal gas. And as expected, when the interaction becomes more and more 
attractive the density profile becomes narrower since the bosons are forced closer together. 
The overall behavior coincides very well with other calculations on condensed bosons [64,74] 
and with experimental determined profiles [56]. This is convincing proof that the correct 
interpretation of the SVM results has been made. 

Correlations 

In accordance with the analysis of the three-boson system, the correlation energy of the 
lowest BEC state in the = 4 case, Ecorr = ^hjj — Etotah is plotted in figure I^TTl for the full, 
the two-body and the mean-field calculations (the Hartree combination is not considered here, 
since it was proven above to be insufficient for treating correlations). The curves from these 
three data sets indicate a clear discrepancy for large negative scattering length, resulting 
in significantly different correlation energies. Focusing first on the mean-field treatment, it 
seems that this description can account for many of the higher-order correlation effects as 
far as a ~ —0.1 bt. This is roughly the same validity region as for the A^ = 3 case and 
also roughly the same as the CP mean-field validity range established in [67] (see A^ = 3 
discussion above). However, the current mean-field calculations are not extensive enough to 
be conclusive in terms of validity estimates. 

^The resulting matrix elements are easily derived from formula (jC.1.3|l . 
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Fig. 5.10: The radial one-body density function defined in ()5.2.3|) for the BEC ground state of a 
system of four trapped bosons, as a function of the distance from the center-of-mass. The different 
curves correspond to SVM calculations for increasingly negative scattering lengths as indicated, 
and have been normalized by / n{r)dr = N. 

Before leaving the mean-field description in this presentation, it is worth noting that in 
the = 4 calculations, the mean-field correlation energy is lower than the full correlation 
energy for a < —0.1 bt, while the opposite is the case for = 3. This might indicate, that 
the pseudopotential approximation adopted in most Hartree-Fock type mean-field theories, 
are destined to include too much correlation energy in the situations where two-body ef- 
fects are dominant, and subsequently too little in the case where higher-order correlations 
are important. However, since almost all GP mean- field theories consider either repulsive 
interactions or only very small negative scattering lengths (in order to satisfy the barrier 
condition \a\N/bt >~ 0.67), this conclusion cannot be confirmed elsewhere. 

Turning the attention now to the results of the two-body and full SVM calculations, 
the answers to some of the interesting questions asked in the introduction of this thesis 
are revealed. First of all, it is evident from the curves in figure 15. IH that three-body and 
higher-order correlations are an integrated part of the dynamics of the four-boson system. 
In addition, these correlations seem to be increasingly important, at least up to a certain 
point, as the interaction between the bosons become more attractive. 

During the further analysis it is convenient to display the difference between the full 
correlation energy, E^^^l, regarded as the "exact" correlation energy, and the two-body cor- 
relation energy, in terms of the former, since this will illustrate the relative importance of 
the higher-order correlations. Such a graph is available in figure 15.121 and indicates that the 
relative effect of higher-order correlations in the four-boson system, saturates when the scat- 
tering length reaches ~ —0.5 bt, i.e. a length similar to half the size of the trap. At this point, 
however, the relative deviation between the two-body and the full treatments correspond to 
almost 35%. In conclusion one may then state, that in the range where — cxo < a < —0.5 bt, 
three- and four-body correlations contribute approximately one third of the correlation en- 
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Fig. 5.11: Correlation energy of the lowest BEC state, Ef,ec,o, in the = 4 boson system, as a 
function of the scattering length for SVM calculations including different levels of correlation. 




Fig. 5.12: The correlation energy difference, Eloor — E^oor normalized by Elorr, as a function of the 
scattering length for the = 4 system of trapped bosons. 

ergy in the = 4 system while two-body correlations contribute the remaining two thirds. 
For a > —0.1 bt, however, the three- and four-body correlations are negligible. 

5.2.3 System with N = 10, and -0.2 6t < a < 

In order to investigate correlations in many-boson systems, it is off course of great interest to 
consider system with more than four particles. The problem is, that the application of the 
SVM for the full correlated description, is only feasible in practice for systems with iV <~ 5 
(because of the symmetrization requirement, see section Fortunately, this is not the 

case for the two-body description where the corresponding computations are independent of 
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Fig. 5.13: Left: Energy levels of the two lowest BEC states for the = 10 system of trapped 
bosons as a function of the scattering length. Right: The corresponding number of lower lying 
(with respect to the BEC state) molecular-type states. All data points refer to the two-body 
correlated SVM calculation. 

N. The next best option in this work is therefore to focus on possible similarities between 
the two-body correlated treatment of the N = 4 system and, for example, the N = 10 
system. 

Energy levels 

The energy levels for the lowest two BEC states of the N = 10 system has been calculated 
in the two-body description for scattering lengths in the range —0.2 6^ < a < (^), as shown 
in figure I^TT^ When looking at the resulting curves, one immediately recognizes the same 
overall behavior as in the three- and four-boson systems. The only clear, but expected, 
discrepancy is in the number of lower lying molecular-type bound states, as can be observed 
on the right of the figure 

5.3 Additional remarks about the results 

The study of boson systems presented in this chapter raises several discussion points, where 
the most vital are concerning the accuracy of the data obtained and, assuming this is accept- 
able, the validity of the interpretation, especially in terms of generalizing the conclusions to 
many-boson systems with > 4. These questions are addressed in this section. 



^The SVM convergence is severely complicated by the fast growing number of lower lying states (as 
clarified in section I^TSjl . which means that the TV = 10 calculations have to be limited to a > —0.2 bt- 
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Fig. 5.14: The energy convergence of the lowest BEC state for the = 4 system of trapped bosons 
in the case where a = —0.13 6f Both the full SVM calculation and the two-body SVM calculation 
is shown. 

Accuracy of the SVM results 

As the observant reader might have noticed, that no comments or indications have been 
made so far about the accuracy of the energies calculated in this work. The reason for 
this is, that the stochastic nature of the SVM makes it inherently difficult to estimate error 
bars on the results since these are all variational upper bounds. In other words, there is 
no definite way to determine how close the stochastic optimization procedure has come to 
reproducing the "exact" wave function. However, one can get a general idea of the accuracy 
achieved in a given calculation from studying the amount of energy gained when adding an 
additional basis function to the basis. 

To exemplify, consider the BEC energy convergence graph displayed in figure l^!T^ where 
the iterative trial and error strategy is illustrated for both the two-body and full calculations 
in the specific N = 4 case where a = —0.13 bt- As a first impression, one notices the several 
dramatic falls followed by abrupt jumps that will occur each time a new lower lying molecular 
bound-state has been found. This behavior is a consequence the specific algoritm employed 
since it automaticly selects a higher target state if the current target state becomes bound. 
Overlooking such "resonances", the overall tendency in the curves is an initially fast and 
subsequently slow decline towards the optimal variational energy that can be achieved with 
the respective trial functions. The key assumption is then, that focusing on the energy gain 
from the last basis increase, and after making sure that this is not within a resonance, one 
can roughly estimate the accuracy of the result from its size. 

Following this example, the accuracy was accordingly checked for the calculations pre- 
sented here. In conclusion to this process, it seems convincing to the writer, that the results 
are accurate to at least three significant digits, which should be enough in the current con- 
text. 
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The molecular-type states 

Although it is the high density molecular-type states that are expected to lead to the 
strongest interparticle relationships, these states are not feasible in a discussion of many- 
body correlations in the current context. The reason is, that the particles in such states are 
so close that they cannot be characterized as weakly interacting and independent of the exact 
shape of the interaction potential, and this means, that the adoption of the Gaussian model 
potential for the two-body interaction, as opposed to one with a hard-core dependence, is 
inappropriate. Furthermore, one can be inclined to think that the influence of correlations in 
the molecular-type systems is so pronounced, that is makes no sense to try and distinguish 
between two-body, three-body etc. effects. In other words, an accurate treatment of the 
molecular-type many-boson systems is inherently tied with the full correlation description. 

Generalization to many-body systems 

In the numerical calculations presented in this chapter, most attention has been given to the 
system of four bosons in a trap. The data obtained for this system is both extensive and 
accurate, and reveal new important information about correlations in a four-body system. 
However, in reference to the introduction of this thesis, it still remains to determine whether 
the conclusions drawn for the four-boson system can be generalized to systems with > 4. 

As is apparent from the analysis of the results, the SVM has limited usability in the 
investigation of large A^ systems. During the restricted time period of this study, the only 
attempt to compare the four-boson results with the features of larger systems, has been 
with the two-body treatment of the A^ = 10 case presented in section 15.2.31 From this effort 
it is clear, that there are definite similarities between the A^ = 4 and the A^ = 10 energy 
levels, at least in the limited range considered. However, the lacking ability to do the same 
comparison in the case of the full correlated treatment, is a severe drawback. 

In conclusion, the stand point taken by this writer is, that the available comparison 
does not constitute enough evidence to allow the interpretations, regarding the effects of 
many-body correlations, to be generalized to A^ > 4 systems as they are. In particular, the 
estimate of the contribution from two-body correlations to the total correlation energy (of 
two thirds) seems inappropriate. It is more reasonable to expect this fraction to fall off as 
the number of particles increase because the number of possible many-body correlations is 
higher. However, the question whether higher-order correlations play an important part in 
the dynamics of many-body systems, has, for all A^, been clearly answered (yes!) from the 
four-body results. 
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Chapter 6 

Conclusion and outlook 



In this thesis the intricate nature of correlations in many-body systems has been studied 
theoretically. The numerical results have been obtained with the Stochastic Variational 
Method using several different forms for the variational trial function. 

First the ab initio theoretical framework of the SVM was derived from the foundation of 
the time-independent Schrodinger equation and an N-body Hamiltonian. It was shown how 
to include different levels of correlation in a variational description by incorporating inter- 
particle dependencies explicitly in the functional form of the trial function. The particular 
descriptions treated were the uncorrelated Hartree-Fock description, the pseudopotential 
mean-field description and the explicitly correlated description. While the first two corre- 
spond to common mean-field theories, the latter includes correlations beyond the mean-field 
and can be designed to introduce only pair-, triplet- or any given higher-order correlation. 
This makes such a description ideal for investigating the significance of correlation effects. 

The SVM was subsequently developed to work with the different correlation descriptions 
by expanding the trial functions in a basis of contracted Gaussians or, for the three-body 
case, in a basis of exponential ftinctions. All the necessary matrix elements were derived in 
the cases where the two-body interaction is given by either a finite Gaussian potential, a zero- 
range pseudopotential or a Goulomb potential. A particular elegant expression was obtained 
for the Jastrow-type two-body correlated trial function, resulting in matrix elements that 
where independent of the number of particles. The trial and error optimization procedure, 
which is the heart of the SVM, was implemented in C++ and thoroughly testet. The 
original few-body algoritm [1] was upgraded to produce fast convergence for systems of 
trapped bosons which require a random value interval over several orders of magnitude. 

After some initial bench marking by extensive calculations on the °°He atom, the main 
part of the numerical work was concentrated on the systems of three and four weakly inter- 
action bosons in a spherical symmetric trap. Attractive interactions corresponding to s-wave 
scattering lengths of — cxd < a < were considered, in the attempt to simulate the conditions 
of experiments on ^^Rb near a Fcshbach resonance. The results obtained reveal the detailed 
behavoir of the energy levels as a function of a, and confirms the gross features found by 
others for the N — 3 system [75]. In addition, both molecular-hke and gas-like states are 
included, whereas standard mean-field methods only treat the latter states. The plotted en- 
ergy curves describe a smooth non-diverging behavior and provides a detailed picture of the 
corresponding physics in combination with an analysis of the effective boson-boson potential. 

Several interesting conclusions were drawn from the numerical results. First of all, it 
was apparent that the mean-field pseudopotential treatment is insufficient for large negative 
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scattering lengths even for few-boson systems. However, the current study was not extensive 
enough to consider the vahdity criteria in detaiL Secondly, the next best treatment corre- 
sponds to assuming that one pair in the system is close in space, but even this two-body 
correlated case eventually fails to reproduce the correct correlation energy, although the 
error is at most 35% in the four-boson system. Subsequently, one would expect this error 
percentidge to climb for system of more particles. In other words, it is not possible to study 
many-body systems accurately over a wide range of scattering lengths, if only two-body 
correlations are taken into account. 

Because all the lower lying molecular-type bound states have to be determined in the 
SVM, the running time of some calculations came close to 12 hours. Therefore, and due to 
the lack of time, it was not possible to do widespread investigations of the = 4 system for 
the full correlated treament. Hence a more systematic investigation of the large negative a 
region could be carried out in a future study, as this scattering range represents a "blind" 
spot in figure I^TT^ which might include additional physics. 

Another immediate extension of the present work is to investigate the five-boson system 
for the two-body and three-body correlated treatments, at least in a limited range of scat- 
tering lengths. Such calculations would represent feasable computations, and the produced 
results can further clarify the relative importance between two-body and higher-order cor- 
relations in an N-body system. Other direct improvements include the adoption of more 
realistic two-body potential models, for example, the simple sum of two Gaussians, where 
one corresponds to a repulsive hard core. For people devoted to computer science, the new 
implemented code also opens the opportunity of rigorous experimentation with the trial and 
error algoritm. This might lead to the discovery of new and better stochastic optimization 
techniques which are becomming increasingly important in computational physics. 

In conclusion, the present investigation of correlations yielded insight into the dynamics 
of many-body systems by investigating the systems of three and four bosons in a trap. The 
results build conclusive evidence of the assumption that higher-order correlations play an 
important role in many-body systems, and that such interparticle mechanisms should not 
be neglected in future studies of this field. 



Appendix A 

Angular momentum functions 



In this appendix, it is shown how to treat systems having non-zero definite angular momen- 
tum (L 7^ 0, 5* 7^ 0) in a variational approach with the general variational trial function 
fj2.6.H) . The goal is to have a the trial wave function that described a system of particles 
having individual orbital angular momenta, spins, Sj, and isospins, tj. The problem is to 
find a set of angular momentum operators that commute with the Hamiltonian and with 
each other, so that common eigenfunctions exist. The corresponding good quantum numbers 
are adapted by the orbital, spin and isospin parts of the basis function. 



A.l Orbital angular moment 

Considering orbital angular momentum first, the single-particle operators, Ij = —irj x V, do 
not commute with the kinetic energy term of the many-body Hamiltonian ()2.1.2|) . However, 
since this Hamiltonian has no relativistic terms, the sum of orbital angular momentum 
operators, L = J2iLi h, does commute with H [7], i.e. 

f L^^ = L(L + 1)^, L = 0,1,2,... 

[H,L] = 0, ^ ^ ' (A.1.1) 

making the definite orbital angular momentum L and projection Ml good quantum numbers. 
The common eigenfunctions of the single-particle operators and kz is the surface spherical 
harmonics, Yi^m^i^'i)- Generalizing, the angular part of the basis function, Olml{,x), is a 
vector-coupled product ^ of spherical harmonics 



^123^^123 



® ■ ■ ■ (^Yi^mA^N)] (A.1.2) 
N 

E Cu\[Yi^^X^;), (A.1.3) 

fe={mi,m2,--- ,mjv} *=1 

where Ck is a product of Clebsch-Gordan coefficients 

Ck = (/lml/2"^2|/l^2-^l2"^l + ni2){Li2mi + m2h'm^\Li2hLi2zmi + 1112 + m^) 

■ ■ ■ {Li2-N-imi + ma ■ ■ ■ mN-ilNfnN\Li2...N-ilNLML) (A. 1.4) 



^Theory for addition of angular momentum can be found in [3]. 
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A. Angular momentum functions 



In this way, each relative orbital motion is being assigned a definite angular momentum and 
GhMLi.^) is dependent on the specific set of angular momenta, {/i, h, ■ ■ ■ , In', L12, L123, ■ ■ ■}, 
chosen. The intermediate momenta, L12, L123, etc., do not in general constitute good quan- 
tum numbers. Thus, for a realistic description, it is often necessary to include several 
different channels, i.e. sets of angular momenta, which is then referred to as the method of 
partial wave expansion [1]. 

The various possible partial wave channels obviously increase the size of the basis. In 
addition, this form of Olmj^ {x), demanding the coupling of (A^— 1) angular momenta, becomes 
increasingly complicated as the number of particles goes up. With a variational approach 
all this can be avoided by adopting another, closely related ^, choice for Oiml{^), proposed 
by Varga and Suzuki in [29]; 

N 

9lMl{x) = \v\'^^YlmA^), with v = ^UiXi = uFx (A.1.5) 

i=l 

Only the total orbital angular momentum enters in this expression. The real vector = 
{ui, . . . ,un}, defining a linear combination of the relative coordinates, may be considered 
a variational parameter, is a positive integer, most often small ^. Like before, several 
terms like ()A.1.5|) in the angular part of the basis function will improve the description. The 
continuity of the parameters Ui yield continuous changes in the evaluated energy expectation 
value, £. This can be more advantageous in a variational calculation than assigning sets of 
discrete angular momenta. 



A. 2 Spin angular momentum: XsMs 

The spin and isospin angular momentum is treated in the same way as the orbital angular 
momentum. If H has no spin terms, the single-particle spin operator, Sj = ^o'i, commutes 
with H. However, only symmetric operators commute with every permutation operator, P, 
used below to ensure the proper symmetry. The symmetric operators 5*^ and Sz correspond- 
ing to the total spin, S = J2iLi'^i, ^^us convenient. The spin part of the basic functions 
then depend on the good quantum numbers S and Ms, and is given by successively coupled 
single-particle spin functions; 



XSMs 



The set of spin quantum numbers, {si, S2, ■ ■ ■ , sn] Su, S123, . . .}, specifies the particular 
coupling. Again, several sets may be needed to obtain a good wave function. Alternatively, 
one can use a spin function based on continues parameters (see [1], sec. 6.4). The isospin 
function, //ta/t' "^^^ be constructed in exactly the same manner as the spin function, Xsms- 



^Any function of the form (|A.1.2|I can be written as a linear combination of terms like HA.1.5|I and vice 
versa. See the detailed proof in [1], sec. 6.2. 

■^Using > in (|A.1.5p corresponds to including many higher partial waves in the expression (|A.1.2|) 
for OlmAx), [29]. 



A. 3 L-S coupling 
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A. 3 L-S coupling 

As already mentioned, the above angular momentum description is valid under the assump- 
tion that the Hamiltonian does not contain any relativistic terms. For the systems under 
consideration in this thesis (see chapter a non-relativistic treatment is sufficient. However, 
this is far from always the case. Accurate calculations of atomic energy levels have to account 
of atomic fine structure, generated by the prominent spin-orbit term, Xlili '^('"i)'^* ' [2]- 
The nucleon-nucleon interaction also has a strong spin-isopin dependence ^. Obviously, nei- 
ther L nor S commutes with the spin-orbit term, while the vector sum, J = L + S, does. 
Consequently, J and Mj will serve as good angular momentum quantum numbers in the 
presence of relativistic terms. The orbital and spin angular momenta are coupled ^ , by 
applying the Clebsch-Gordan series to the L and S quantum numbers of the separate parts, 
&lMl{x) and Xsms: &s indicated in the form of the basis functions ()2.6.1|) . 



''Modern two-body NN potentials are Argonne vis, Nijmegen II, Reid93, CN-Bonn [8]. 
^Known as Russell-Saunders or L-5 coupling, [7]. 
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Appendix B 

Hartree-Fock ground state of identical 
fermions 

As Wolfgang Pauli's famous exclusion principle states, identical fermions cannot occupy the 
same quantum state at the same time. This means that in an idealized case at T = the 
Hartree-Fock many-fermion ground state, ^E^^^, should be formed by occupying N single- 
particle energy levels from the lowest up. If all interactions where neglected the many-fermion 
state would then represent a filled Fermi sphere where all stationary orbitals corresponding 
to an energy less than the Fermi energy [Ep = ksTp) is occupied by exactly one particle. 
However, interactions pertube this picture and modify the energy levels and single-particle 
states. 

Within the Hartree-Fock method this modification is well described by the variational 
single-particle wave functions. Assuming that the Pauli principle is still forced on the product 
wave function 1)3.2.11) by making it explicitly antisymmetric ^ and that the N single-particle 
states 0102 • • • 0Ar with lowest energy are given by the spin-orbitals 

0,(r) = ^.(r)xi/2,s„ z = 1, 2, . . . , iV (B.0.1) 

where Si is the spin of the zth fermion with < Xi/2,si\xi/2,sj >= ^SiSj, the interpretation of 
the Hartree-Fock equations as eigenvalue problems [2] with the Lagrangian multipliers. Si, 
as the one-particle eigenvalues, yields solutions 

N 

Si =< ipi\hi\ilJi > + i^ii^jlVijlipiipj - Ss.s.'ipjipi > (B.0.2) 
i=i 

These eigenvalues obey Koopman's theorem [57], Ehf{N) — Ehf{N — 1) = Sn- Then Si is 
the energy needed to remove the ith fermion from the system provided the change in the 
wave function for the other particles can be neglected (e.g. when N ^ 1). Summing over 
Si, and comparing with eqs. ()3.2.2|) - ()3.2.4|) . the total energy of the fermion ground state in 
the Hartree-Fock approximation is not J2iLi ^-s might be expected but rather 

TV N 

EPf =< ^SfI^I^Sf >=Y.^^-Y.< ^^^^V^^^^^J - Ss^s.^A > (B.0.3) 

i=l i<j 

^E.g. by adopting the Slater determinant description, ^/fi? — A'i'H — det \(j)i<j)2 ■ ■ ■ (pw], [7]- 
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B. Hartree-Fock ground state of identical fermions 



The "additional" term can be understood as eliminating the double counting of pairs of 
particles since £i includes the energy for each particle interacting with every other 

particle and so counts the contribution from a given pair twice. 

Equation ()B.0.2|) illustrates the influence of fermion-fermion interactions on the single- 
particle energy levels. This effect is well described in the Hartree-Fock approach. However, 
another effect of interactions is that the fermions might be scattered in and out of the single- 
fermion levels, which are no longer stationary. This is not supported within the independent 
particle approximation. Fortunately, a fermionic particle can only change energy in a collision 
if the final energy state is unoccupied. In the ultra-cold quantum regime, it is highly likely 
that low-energy states are already occupied. Thus the validity of the Hartree-Fock method 
for the fermionic ground-state relies on Pauli blocking since it tends to suppress any process 
in which fermions change energy states. 



Appendix C 
Matrix elements 



C.l Explicitly correlated Gaussian basis 

All matrix elements of the correlated Gaussian basis functions ()4.2.3p are given by integrals 
that can be evaluated to simple analytical results. Since the systems considered in this 
thesis are limited to central interactions and zero angular momentum, all that is needed to 
calculate the matrix elements are the three basic integral formulas [1] 

3 

poo poo //o_\7V— 1\ 2 



oo J ~oo 
oo /•oo 



/oo poo 

■■■ dxidx2 ■ ■ ■ dxN-iX^Cxe-^^'^'^^'^ = 3Tr(A-^C)Jo (C.1.2) 
■oo J ~oo 

/oo poo , „ ^ 1 

■J dxidx2---dxN-iS{c^x-r)e-^^''''^'' = i^—ye-^''-'lo (C.1.3) 



where x'^ = {xi, X2, . . . , x^_i) are independent coordinates, A and C are symmetric positive 
definite matrices and c^^ = c^A^^c. 

Introducing B = A*-^ ^ + A*-^-* for notational convenience and using transformations 
ri = [u^'^'^Y'x and Vij = {u^'^^^)'^x, defined in (I4.2.4j) . the above integrals give the following 
analytical expressions for the overlap, kinetic, trap and interaction matrix elements respec- 
tively 



1 - 



1 



--Tr(S-U(^'')AAW)(^,#,) 



(C.1.5) 
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^1 3 ^ 



i=l 



i=l 
N 



i=l 



(C.1.6) 



/oo 

where in the last expression the dependency on the specific form of the central potential 
describing the interaction, V{r), is isolated in the single integral 



„ . 3 ^OO 



t'(%) = (g:j^ / rfrl^(r)e-4^-^^ (C.1.8) 



oo 



In the case of equal masses, rrii = m, one may simplify the trap expression ()C.1.6|) with the 
identities 

N 

= I and it^*^) = 2 (C.1.9) 

i=l 

that are satisfied by the Jacobi transformation. To complete the evaluation of the matrix 
elements the v{cij) expressions for the three simple central potentials used in this work are 
listed: 

• For few body atomic systems the Coulomb interaction, V{r) = is used where qi 
and q2 are the charges. In this case the integral ()C.1.8|) gives [18] 

vcouiombicij) = ^irqiqj (^^) ' ^ dr re"^^"^'-" = 2qiqj,j^ (C.1.10) 

• In the calculations of N-body boson systems the interaction is described by a Gaussian 
potential, V{r) = Vqc^^ , giving 

s 3 roo 1 3 

.c..„(e.,) ^ V,(|)'/^*e-J<--^'- ^ V.(^:^)' (C.l.ll) 

• The commonly used "mean-field" description of Bose- Einstein Condensates has a two- 
body interaction given by a zero-range delta function potential [52], V{r) = ^5(r), 
in which case 

VDeltaiCij) = (C.1.12) 



where a is the s-wave scattering length. 
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C.2 Two-body correlated Gaussian basis 

In this section, the matrix elements {i'k'\Hint\4'k) and (i'k'li'k) for the two-body correlated 
version of the explicitly correlated Gaussian basis functions defined in ()4.5.1()|) . are evaluated. 
Introducing 'y^'^^ = j3'^^'> — a^^'^ for notational convenience, the basis functions written in Jacobi 
coordinates are given by 



1 ^ 



k ' 



4^' 



where 



{C.2.1] 



(C.2.2) 



Since we are only considering systems of identical particles the Hamiltonian, Hint, is sym- 
metric, and one can initially apply ()4.3.5|) to obtain 



N 



i<j 



and likewise for the overlap matrix element {Hint !)• Using the previously derived matrix 
element formulas ■4jl - ()( ^1 .7;) with rrii = m, the expressions for the overlap, kinetic, trap 
and interaction matrix elements between V'i'^'' and are simply 



(2^) 



N-l \ f 



detS 



Af-l 



i=l 
^ 1 



3fi 



i=l 



N 



-mu) 



N 



'■1x(B-'){i.f'\i.f^) 



1 



(C.2.4) 
(C.2.5) 

(C.2.6) 

(C.2.7) 



where B = A*^^^'*"' -\- A^^-''^\ Such (naive) adoption of the general formulas for the evaluation 
of the two-body matrix elements apparently leads to N different contributions to ()C.2.3|) for 
the overlap, trap and kinetic terms and even N'^{N — l)/2 for the interaction term. This 
is however not the case since most of the terms are identical. Consider the explicit Jacobi 
coordinate representation of i^f^^ and defined by the matrices 



^(12,fc) 



Na^'''> 











\ 




(C.2.8) 
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^(13,fc) 







f 7('=) + A^aW 
iVa('=) 



v 















^(23,fc) 







|7('=) + Na^''^ 
iVa^^) 



\ 




\ 






(C.2.9) 



(C.2.10) 



and 



^(34,fc) 







V 





7^'=) + A^a^'^) 
('^^ + A^a^'^) 








y2|7(^) + ATaC^) 














\ 




A^aC^) I 



(C.2.11) 



Using these matrices and the general definition ()C.2.2|) one may convince oneself ^ that both 
detS,Tr(S-^), Tr(S"^ A^^^''^'^ A^^^''^)) and even u^^""^^ B-^u^'^''\ will evaluate to only three 
different values and that this, in turn, allows the very important simplification 



t = l, 3 = 2 
1 = 1,2, J = 3, 
^ = 3,...,A^, J 



.,N 

4 AT 



(C.2.12) 



and correspondingly for the overlap matrix elements. Moreover, the A^(A^ — l)/2 terms in 
the sum of the interaction matrix element are also limited to a constant number of different 
values. The evaluation of these terms is strait forward but rather extensive and here only 
the end results will be listed: 



N 



{12)\ 



v{cu) + 2{N - 2)v{cis) + ^{N -2){N- 3)v{cu) 



(C.2.13) 



N 



i'^i \ Yl ^n^n\rk) = [^(^12) + ^^(Cis) + t^(c23) + (A^ - 3){v{Cu) + ^^(034) + ^^(034)) 



1 



+ -(Ar-3)(Ar-4Mc45) x (V^f )|V^f )) 



(C.2.14) 



^E.g. using [9]: 



a b 
b c 



= ac — bb, and 



a b 
b c 



-bb 



C.3 Correlated exponential basis (N = 3 only) 
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N 



(V-f^l Vmnl^pt^) = ^(Cl2) + MCls) + v{c2a) + 1^(034) + 2{N - 4) (i;(ci5) + ^^(Css)) 



+ -{N - 4){N - 5)v{c,,)\ X (^f^l^f)) 



(C.2.15) 



where it is assumed that > 4 (or that only the appropriate terms that would be available 
if < 4 are taken into account) and = u^'^'^)'^ B~'^u^'^'^\ The interaction potential 
function, f (cnm), has been derived in ()C.1.10|1 - ()C.1.12|1 for the central interactions of interest 
in this thesis. 

Collecting the above results, the matrix elements of Hint for the two-body correlated 
basis functions (and similarly for the overlap elements using Hint 1), can be written as a 
sum of three terms, i.e. 



(V'fc'l^mil^fc) 



2 

m 



(^f + N,,{i;lY>\Hint\i^r) + N,,{i;lY>\H,nt\i^D 



',(12)\ 



.,(12) 



(13) N 



',(12) I 



',(34)\ 



(C.2.16) 

where A^is = 2(A^ - 2) and A^34 = (A^(A^ - l)/2 - 1 - 2{N - 2)) and, more importantly, each 
of the individual terms is given by a combination of a few expressions with computational 
complexity 0{1). 



C.3 Correlated exponential basis {N = 3 only) 

In this section, the explicit analytical expressions for matrix elements needed in the exponen- 
tial basis variational solution of a nonrelativistic Coulombic three-body system with L = 0, 
are presented. Additional formulas for arbitrary values of total angular momentum can be 
found in [44]. 

Evaluation of the elements of the overlap matrix (j2.4.7|) and the Hamiltonian matrix 
()2.4.6p corresponding to the Coulombic three-body Hamiltonian (in a.u. units) 

= -iv^AV. + «1* + + «L ,c.3.1) 

2 Xi X2 Xs 

where x'^ = (ri2, Vy^, r23) and A is defined by (j2.3.9j) . requires only the elements 

(V'fc'l^^fe), {M-\§i^k), {i^k'\^a.-^a,,\SiJk), (C.3.2) 

Xi 



to be calculated for all z, j = 1, 2, 3 and 5 = ^(1 + P12 + P13 + P23 + ^12^^13 + ^^12^^23)- With 
the simple exponential basis function, ipk = exp(— a^xi — PkX2 — 'Jk^s), the scalar integral is 
advantageously defined by [46] 



( ) = / / / ri2ri3r23dri2dri3dr23 = XiX2X3dxidx2dx3 (C.3. 3) 





even though the interparticle distances are not independent variables. In this case, however, 
the elements in fjC.3.2|) are conveniently expressible in terms of the basic three-body integral 

F{ni,n2,n3) = J J J Xi^X2^x^^dxidx2dx3 exp{—axi — Px2 — 'JX3) (C.3. 4) 
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where a = ak' + ak,P = Pk' + Pk and 7 = 7fc' + 7fc. Simple analytical formulas for all necessary 
F{ni, n2, ris) integrals are derived at the end of this section. 
The overlap integral in the above notation is then simply 

(^fc'l^fc) ^ J J J xi^2X3dxidx2dx3 exp{-axi - (3x2 - 7^3) 

= F(1,1,1) (C.3.5) 
Correspondingly, the matrix elements of the potential energy terms become 

(C.3.6) 
(C.3.7) 
(C.3.8) 

To determine the matrix elements of the kinetic energy terms the gradient operator, V, is 
needed. In the relative coordinates it has the form [44] 

The angular part in the gradient does not contribute when working on purely radial basis 
functions, hence 

xf + x'j-xl 92 

v.. V.^. = <j '"^'^^^ ^ '"^'''''^ (C.3.10) 



{A' 


1 

Xi 


m 


= F(0,1,1) 




1 

X2 




= F(1,0,1) 




1 

X3 




= F(1,1,0) 



dx^ Xi dxi ' 

where the law of cosines [18], Xi- Xj = ' , , was used. The kinetic energy elements are 
then 

(V'ik'l V' \^k) = {^^2 + -^\^k) = -alF{l, 1, 1) + 2akF{Q, 1, 1) (C.3.11) 



asi I 

052 I 



dxl 


+ 


Xx 


dxi 


Q2 


+ 


2 


d 


dxl 


X2 


dX2 
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+ 


2 


d 



iV^fc) = -PiF{l, 1, 1) + 2/5fcF(l, 0, 1) (C.3.12) 



{M'^IM = {M^ + -^\^k) = -llF{l, 1, 1) + 27feF(l, 1, 0) (C.3.13) 



[F(2, 0, 1) + F(0, 2, 1) - F(0, 0, 3)] (C.3.14) 



-]^akik [F{2, 1, 0) + F(0, 1, 2) - F(0, 3, 0)] (C.3.15) 



C.3 Correlated exponential basis (N = 3 only) 
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^ _|_ ^2 ^2 

= -iA7i[-F(l,2,0) + F(l,0,2)-f(3,0,0)] (C.3.16) 

Since the permutation of particles does not change the exphcit form of ipk (e.g. Pui^k = 
exp(— — akX2 — Jk^s)), one can determine all terms in the matrix element expressions 
m (frT3:2|l from the ab ove results. 

Analytical formulas for the basic three-body integral 

To complete the derivation of the three-body matrix elements a few cases of the basic 
three-body integral F{ni,n2-,n^) in ()C.3.4|) has to be calculated. Transforming to truly 
independent perimetric coordinates given by 

Ui = ]^{rik + rij-rjk), i ^ j ^ k = {1,2,?,), (C.3. 17) 

the integral F(0, 0, 0) is trivial [46] 

Other cases of F(?t,i, ?t,2, ^3) can be derived by differentiating or integrating this expression 
with respect to a, 13 and 7. Introducing the function 

D{a, b, c) = ^, ,^ ^ , (C.3.19) 

the necessary F{ni, ^2, jt-s) integrals are given by 

F(l, 1, 0) = D{2, 2, 1) + D{2, 1, 2) + D{1, 2,2) + 2* D{3, 1, 1) (C.3.20) 

F{3, 0, 0) = 2 * 3 * [D{3, 1, 2) + D{2, 1, 3) + D(4, 1, 1) + D{1, 1, 4)] (C.3.21) 
F(l, 1, 1) = 2 * [D{3, 2, 1) + D{2, 3, 1) + D{1, 3, 2) + D{1, 2, 3) 

+ D{2, 1, 3) + D{3, 1,2) + D{2, 2, 2)] (C.3.22) 
F{2, 1, 0) = 2 * [D{3, 2, 1) + D{2, 1, 3) + D{1, 2, 3) + D{2, 2, 2) 

+ 2*L)(3,1,2) + 3*L>(4,1,1)] (C.3.23) 

Permutation of the coordinates Xi X2 ^ X3 in the basic integral ()C.3.4|) corresponds to 
permutation of {ni,a) ^ {n2,(3) ^ {n2,'y) and allows easy construction of the remaining 
cases, e.g. F(0, 1, 2) = F{2, 1, 0; a ^ 7). 
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Appendix D 



CH — h implementation of the 
Stochastic Variational Method 



As with most numerical calculations in physics the main effort during implementation is 
on precision and speed. With this particular method, heavy duty demands during matrix 
element calculations and large eigenvalue problems, makes an efficient routine essential. In 
this section, considerations on key aspects of implementing SVM are described. 

D.l Arbitrary precision arithmetic 

When working with a very large basis the standard 64-bit precision arithmetic might not be 
sufficient to maintain numerical stability in the computations. To this end, the free multi- 
precision package doubledouble [41] is applied for calculations with ffoating point numbers 
of an 128-bit length. A C++ class was then wrapped aroud this type to allow for 64-bit 
exponents. However, employing this class makes all computations aproximately 7 times 
slower. 

D.2 Scaling overlap values to minimize loss of accuracy 

The magnitude of the overlaps is the dominant scale of the matrix elements corresponding 
to a given trial function. Unfortunately, scaling the overlap to (ipk'li^k) ~ 1 is not possi- 
ble without breaking up the (binary) power function calculation and scaling concurrently. 
However, one can use the overlap magnitude to estimate the maximum and minimum values 
handled in the eigenvalue solution and hence set up a validity check for a calculation on a 
specific computer. Neglecting all constant factors ^ in the overlap expression ()C.1.4|) while 
introducing a scale factor S gives 



^For a running calculation TV can be viewed as a constant factor, although, when no scaling is applied it 
must be kept in the power expression to avoid overflows when computing the {N — 2)th power for very large 
N (i.e. TVa ~ 1). 




rnm 




(D.2.2) 
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D. C++ implementation of the Stochastic Variational Method 



where amax and amin is the maximum and minimum value possible for a and correspondingly 
for 7. One can center the overall magnitudes around 10° by multiplying all overlaps with 

the factor S = ^ O^maxCy-min- 

D.3 Avoiding linear independence 

At the core of SVM is the random trial and error selection of the basis functions \E'(aj),i = 
1..K. The result is a state-space spanned by a finite number of dense nonorthogonal func- 
tions. Because of the random origin of the basis functions they cannot be expected to be 
linearly independent. Although linear dependence is seldom with the fully random basis, in 
contrast to geometric progression and random tempering [1], care must be taking to avoid 
it. In practical problems exact linear dependence, like degeneracy in the basis, is unlikely, 
but still close to exact linear dependence between basis functions will lead to poor precision 
in the calculation. This is because one or several eigenvalues of S gets very small when the 
linear dependence is distinct, producing large expansion coefficients in the trial function. 
Then a small error in calculation of the matrix elements of H and S can result in a large 
error in the variational energy. 

D.4 A symmetric-definite generalized eigenvalue prob- 
lem 

Adding a new trial function to the basis demands solving a symmetric-definite generalized 
eigenvalue problem with good precision and efficiency. Solving eigenvalue problems has 
been an intense area of research since the dawn of computers in the 1950's, resulting in 
numerous elegant methods, specially designed for different conditions of the eigenproblem 
(see summary of research in [17]). Fortunately, for real symmetric-definite matrices, the 
eigenproblem is relatively simple. The eigenvalues are always real and there exists a complete 
orthogonal eigensystem that is exploited in very efficient numerical methods. 

For real symmetric matices, the eigenproblem is relatively simple, due to the existence of a 
complete orthogonal eigensystem, and the fact that all eigenvalues are real. These properties 
are exploited in the most efficient numerical methods, and the symmetric eigenproblem may 
be considered as solved: for small matrices n <= 25 we have the QR method, one of the most 
elegant numerical techniques produced in the field of numerical analysis; for large matrices 
25 < n < 1000, we have a combination of divide and conquer with QR techniques. For 
asymmetric matrices the picture is less rosy. 

D.5 Root finding 

The Gram-Schmidt orthogonalization formula (j4.1.2|) implies that 




(D.5.1) 



D. 6 Making sure A is positive definite 
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with 

K 

Uk+iW = {^{aK+i)maK+i)) - Yl I(*(«i^+i)l0.>r = (MaK+iMK+i) (D.5.2) 

i=l 

and 

K+l 



i=l 



i=l 



j=l \ i=l 

K 



i=l 



Hence = {(l)j\H\(j)K+i) for j — 1..K and the most efficient way to implement the 

root finding is by using the expression 



1=1 ^ ' 



= V ^ ^ (^if+il^il^i^+i) - A||0i^+i|| 

1^1 ~ ^) 

Notice that this equation also proves that the energy will be lower when the dimension of 
the basis increases. 



D.6 Making sure A is positive definite 

An N X N real matrix A is called positive definite if 

x^Ax > 0, 

for all nonzero vectors x e R^. There are various ways to test if a matrix A is positive 
definite based on the following observations [16]: (a) all the eigenvalues of a positive definite 
matrix are positive, (b) all upper left (i.e. principal) submatrix determinants are positive 
and (c) a real symmetric matrix is positive definite iff there exists a real non-singular lower 
triangular matrix L such that 

A = LL^. (D.6.1) 

The latter approach, called Cholesky factorization, is the most efficient in the case of a 
large size symmetric matrix. This method is implemented based on the ALGOL procedure 
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Choldetl from [13]. The factorization algoritm fails if the matrix is not positive definite. 
In addition, observation (b) above is used to include simple and fast special cases for = 
1,2,3. For such low A^, only a few simple determinants have to be evaluated, making the 
factorization routine cumbersome in comparison. This will speed up calculations involving 
four or less particles. 

D.7 Inversion of positive definite symmetric matrices 

Inversion of positive definite symmetric matrices is effectively done by the Gauss- Jordan 
Method following the lines of the ALGOL procedure gjdef 2 from [14]. The lower triangular 
array representation of symmetric matrices allow fast element access. 

D.8 Symmetry: all possible permutations 

To implement the symmetrization procedure described in section 14.31 one need to find all 
possible permutations of the set of identical parameters. The SEPA algorithm [53] is used to 
create all permutations from which the corresponding linear transformations of the Jacobi 
coordinates for each trial encountered is evaluated. This is done only once for each trial and 
then stored for later uses. 



Appendix E 

Program usage information 



Usage: scatlen [OPTIONS] 

Calculates the scattering length for a two-body interaction of identical bosons 



-h — help 
-mass <float> 
-pot <int> 
-VO <float> 
-b <float> 
-rmeix <float> 
-steps <int> 
-digits <int> 
-printpot 
-printwave 
-compare 
-notxt 



Prints this usage message. 

Mass of the particle in atomic mass units (default 86.9091835) 
Type of potential =' gauss' or 'square' (default is 'gauss') 
The potential amplitude in a.u. (default is -5.986e-8) 
The potential width a.u. (default is 18.9 = 1 nm) 
Integrate to this radius in a.u. (default 4b) 

Force a specific number of integration steps (default is 10000) 
Force a specific number of correct digits (default 4) 
Prints the potential points as 'rl pi | r2 p2 | . . . ' 
Prints the radial wave function as 'rl wl I r2 w2 | . . . ' 
Compare result with analytical square box value or Born approx 
Demand scattering length as only output 
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Usage: bee [OPTIONS] 

Calculates the energy for a given state of an N-body system using the Stochastic Variat 
ional Method. 



-h — help 
-N <int> 
-state <int> 
-basis <type> 
-par <int> 
-sym [int] 
-antisym [int] 
-size <int> 
-trials <int> 
-times <int> 
-reps <int> 
-recycle <int> 
-rfine <float> 
-rtimes <int> 
-rbegin <int> 
-rend <int> 
-units <type> 
-notrap 
-dotrap 
-int <type> 
-b <float> 
-VO <float> 
-as <float> 
-aB <float> 
-seed <int> 
-rint <float> 
-rmin <float> 
-rmax <float> 
-rbint <float> 
-rbmin <float> 
-rbmax <float> 
-rlog [<int>] 
-fin <name> 
-fout <naine> 
-digits <int> 
-noimps <int> 
-Idep <float> 
-threads <int> 
-save <int> 
-endtime <int> 
-check 
-warn 
-Stat 
-no info 
-notxt 



Prints this usage message . 

Number of particles (default is 3) 

Specify <int>th, 'pos ' , 'neg' eigenstate as target (default is 0) 

Basis type = 'full' , 'twobody' , 'hartree' (default is 'twobodyO 

Override the number of nonlinear parameters in the full basis 

Symmetrize trials [only first <int> particles] (default is N) 

Ant i symmetrize trials [only first <int> particles] (no default) 

Size of the basis to be calculated (default is 10) 

Number of trials pr. nonlinear parameter (default is 100) 

Number of times to restart trials loop at par. one (default is 10) 

Number of times to repeat the trial&error procedure (default is 1) 

Recycle intermission every <int>'th new basis (default is size) 

Recycle type: = random, <float> = finetune (default is 0) 

Number of times to repeat recycle procedure (default is 0) 

Recycle procedure should begin at basis <int> (default is 1) 

Recycle procedure should end at basis <int> (default is K) 

Calculation units = 'hou' or 'au' (default is h.o.u.) 

Remove trap from system (to calculate bound states) 

Add trap to the system (to cancel previous -notrap) 

Interaction = 'non' , 'zero' , 'gauss' , 'coulomb' (default is non) 

Set the potential range in a.u. (default is 11.65) 

Set the potential amplitude in a.u. (default is 1.103130e-7) 

Specify the scattering length in a.u. (used only for output) 

Override the calculated Born scattering length in a.u. 

Seed for the random number generator (default is 1) 

Random interval range for alpha coefficient (default is 10.0) 

Override the estimated alpha random interval minimum 

Override the estimated alpha random interval maximum 

Random interval range for beta coefficient (default is 10.0) 

Override the estimated beta random interval minimum 

Override the estimated beta random interval meiximum 

Use logarithmic random interval [with base <int>] 

Filename for basis input (default none) 

Filename for basis output (default none) 

Number of digits used in rootfinding and output (default is 8) 
Succeeding 'no improving trial's allowed (default is 5) 
Specify the lowest linear dependency allowed (default is le-6) 
Maximum nximber of cpu threads used (default is 2) 
Save basis every <int>'th minute (default is 10) 
Time limit for the calculation in minutes (default is no limit) 
Check explicitly for numerical instabilities (default is off) 
All warnings are displayed (default is off) 
Post-calculation statistics are displayed (default is off) 
Output only calculation results (default is with info) 
Output only: [basisnumber energy] (for use with e.g. MATLAB) 
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-eigenvalues Output all eigenvalues at the end (i.e. show excited states) 
-result Output only the final energy result (for use with e.g. MATLAB) 

-resultall Output for MATLAB by writing information and end results like: 

[N energy as b VO aB aho eho K state Rrms mean (alpha ' s) mean(beta's) dE seconds] 



-bee For Bose-Einstein Condensate calculation (default is Rb87) 

-mass <float> Mass of BEC boson in a.m.u. (default is m(Rb87) =86. 9091835) 
-freq <float> Specify the trap frequency in Hz (default is 77.87) 



-bound For N-body bound state calculation (default is Helium atom) 

-masses <list> Set particle masses ml,m2,..mN in a.u. (default is l,l,le300) 
-charges <list>Set particle charges ql,q2,..qN in a.u. (default is -1,-1,2) 
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